\documentclass[10pt]{article}
\usepackage[verbose, a4paper, hmargin=2.5cm, vmargin=2.5cm]{geometry}

\usepackage{fontspec}
\usepackage{ctex}
\usepackage{paratype}


\usepackage{amsmath}
\usepackage{amsfonts}
\usepackage{amssymb}
\usepackage{esint}

\usepackage{graphicx}
\usepackage[export]{adjustbox}
\usepackage{mdframed}
\usepackage{booktabs,array,multirow}
\usepackage{adjustbox}
\usepackage{tabularx}
\usepackage{hyperref}
\hypersetup{colorlinks=true, linkcolor=blue, filecolor=magenta, urlcolor=cyan,}
\urlstyle{same}
\usepackage[most]{tcolorbox}
\definecolor{mygray}{RGB}{240,240,240}
\tcbset{
  colback=mygray,
  boxrule=0pt,
}
\graphicspath{ {./images/} }
\newcommand{\HRule}{\begin{center}\rule{0.9\linewidth}{0.2mm}\end{center}}
\newcommand{\customfootnote}[1]{
  \let\thefootnote\relax\footnotetext{#1}
}

\begin{document}
\section*{第三章 随机模拟}

\section*{3.1 概述}

在用数学模型, 包括概率统计模型处理实际应用中的问题时, 我们希望建立的模型能够尽可能地符合实际情况。但是, 实际情况是错综复杂的, 如果一味地要求模型与实际完全相符, 会导致模型过于复杂, 以至于不能进行严格理论分析, 结果导致模型不能使用。所以, 实际建模时会忽略许多细节, 增加一些可能很难验证的理论假设, 使得模型比较简单, 可以用数学理论进行分析研究。

这样, 简化的模型就可以与实际情况有较大的差距, 即使我们对模型进行了完美的理论分析, 也不能保证分析结果是可信的。这一困难可以用随机模拟的方法解决。

模拟是指把某一现实的或抽象的系统的某种特征或部分状态, 用另一系统 (称为模拟模型)来代替或模拟。为了解决某问题, 把它变成一个概率模型的求解问题, 然后产生符合模型的大量随机数, 对产生的随机数进行分析从而求解问题, 这种方法叫做随机模拟方法, 又称为蒙特卡洛 (Monte Carlo) 方法。

例如, 一个交通路口需要找到一种最优的控制红绿灯信号的办法, 使得通过路口的汽车耽搁的平均时间最短,而行人等候过路的时间不超过某一给定的心理极限值。十字路口的信号共有四个方向, 每个方向又分直行、左转、右转。因为汽车和行人的到来是随机的, 我们要用随机过程来描述四个方向的汽车到来和路口的行人到来过程。理论建模分析很难解决这个最优化问题。但是, 我们可以采集汽车和行人到来的频率, 用随机模拟方法模拟汽车和行人到来的过程, 并模拟各种控制方案, 记录不同方案造成的等待时间, 通过统计比较找出最优的控制方案。

随机模拟中的随机性可能来自模型本身的随机变量, 比如上面描述的汽车和行人到来, 也可能是把非随机的问题转换为概率模型的特征量估计问题从而用随机模拟方法解决。

例 3.1.1. 为了计算圆周率 \(\pi\) 的近似值可以使用随机模拟方法。如果向正方形 \(D = \{ \left( {x,y}\right)\) : \(x \in  \left\lbrack  {-1,1}\right\rbrack  ,y \in  \left\lbrack  {-1,1}\right\rbrack\) 内随机等可能投点,落入单位圆 \(C = \left\{  {\left( {x,y}\right)  : {x}^{2} + {y}^{2} \leq  1}\right\}\) 的概率为面积之比 \(p = \frac{\pi }{4}\) 。如果独立重复地投了 \(N\) 个点,落入 \(C\) 中的点的个数 \(\xi\) 的平均值为 \({E\xi } = {pN}\) , 由概率的频率解释,

\[
\frac{\xi }{N} \approx  \frac{\pi }{4},\;\pi  \approx  \widehat{\pi } = \frac{4\xi }{N}
\]

可以这样给出 \(\pi\) 的近似值。

随机模拟方法会引入所谓随机模拟误差: 上例中估计的 \(\widehat{\pi }\) 实际是随机的,如果再次独立重复投 \(N\) 个点,得到的 \(\widehat{\pi }\) 和上一次结果会有不同。这是随机模拟方法的特点,即结果具有随机性。因为结果的随机性导致的误差叫做随机模拟误差。

使用随机模拟方法, 我们必须了解随机模拟误差的大小, 这样我们才能够设计合适的重复试验次数来控制随机模拟误差。比如,这个例子中 \(\xi\) 服从 \(\mathrm{B}\left( {N,\frac{\pi }{4}}\right)\) 分布,有

\[
\operatorname{Var}\left( \widehat{\pi }\right)  = \frac{\pi \left( {4 - \pi }\right) }{16N},
\]

由中心极限定理, \(\widehat{\pi }\) 近似服从 \(\mathrm{N}\left( {\pi ,\frac{\pi \left( {4 - \pi }\right) }{16N}}\right)\) 分布,所以随机模拟误差的幅度大约在 \(\pm  2\sqrt{\frac{\pi \left( {4 - \pi }\right) }{16N}}\) (随机模拟误差 95\% 以上落入此区间)。

这样的随机误差幅度也是随机模拟误差的典型情况。一般地,假设随机变量 \(X\) 期望为 \(\theta\) ,方差为 \({\sigma }^{2}\) 。产生随机变量 \(X\) 的 \(N\) 个独立同分布随机数 \({X}_{i},i = 1,2,\ldots ,N\) ,用样本平均值 \({\bar{X}}_{N} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{X}_{i}\) 估计 \(\theta\) ,由中心极限定理可知 \({\bar{X}}_{N}\) 近似服从 \(\mathrm{N}\left( {\theta ,{\sigma }^{2}/N}\right)\) ,于是随机模拟误差幅度为 \({O}_{p}\left( \frac{\sigma }{\sqrt{N}}\right) \left( {N \rightarrow  \infty }\right)\) 。为了估计 \({\bar{X}}_{N}\) 的渐近标准差 \(\sigma /\sqrt{N}\) ,可以用样本方差 \({S}_{N}^{2} = \frac{1}{N - 1}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {X}_{i} - {\bar{X}}_{N}\right) }^{2}\) 代替 \({\sigma }^{2}\) 进行计算。为了控制估计的标准差小于 \({\sigma }_{0}\) ,可以先取较小的 \({N}_{0}\) 抽取 \({N}_{0}\) 个样本值计算出 \({S}_{{N}_{0}}^{2}\) ,用 \({S}_{{N}_{0}}^{2}\) 估计 \({\sigma }^{2}\) ,然后求需要的 \(N\) 的大小:

\[
\frac{{S}_{{N}_{0}}}{\sqrt{N}} < {\sigma }_{0},\;N > \frac{{S}_{{N}_{0}}^{2}}{{\sigma }_{0}^{2}}.
\]

用 \({\bar{X}}_{N}\) 估计 \(\theta  = {EX}\) 时,也可以利用中心极限定理计算 \(\theta\) 的近似 \({95}\%\) 置信区间:

\[
{\bar{X}}_{N} \pm  2{S}_{N}/\sqrt{N}. \tag{3.1}
\]

随机模拟方法虽然避免了复杂的理论分析, 但是其结果具有随机性, 精度很难提高: 为了增加一位小数点的精度,即误差减小到原来的 \(\frac{1}{10}\) ,重复试验次数需要增加到原来的 100 倍。 随机模拟方法有如下特点:

\begin{itemize}
\item 应用面广,适应性强。只要问题能够清楚地描述出来, 就可以编写模拟程序产生大量数据, 通过分析模拟数据解决问题。
\end{itemize}

\begin{itemize}
\item 算法简单, 容易改变条件, 但计算量大。
\end{itemize}

\begin{itemize}
\item 结果具有随机性,精度较低 (一般为 \(O\left( \frac{1}{\sqrt{N}}\right)\) 级)。
\end{itemize}

\begin{itemize}
\item 模拟结果的收敛服从概率论规律。
\end{itemize}

\begin{itemize}
\item 对维数增加不敏感。在计算定积分时, 如果使用传统数值算法, 维数增加会造成计算时间指数增加,但是如果使用随机模拟方法计算,则维数增加仅仅造成不多的影响。
\end{itemize}

随机模拟用在科学研究中, 常常作为探索性试验来使用。假设科学家有了一个新的模型或技术的想法, 但是不知道它的效果怎样所以还没有对其进行深入的理论分析, 就可以用随机模拟方法大量地重复生成模拟数据, 根据多次重复的总体效果来判断这种模型或技术的性能。如果模拟获得了好的结果, 再进行深入理论分析对模型进行完善; 如果模拟发现了这个模型的缺点, 可以进行有针对性的修改, 或者考虑转而其它解决办法。

随机模拟在科学研究中的另一种作用是说明新的模型或技术的有效性。在公开发表的统计学论文中, 已经有一半以上的文章包括随机模拟结果 (也叫数值结果), 用来辅助说明自己提出的模型或方法的有效性。有时因为对模型或方法很难进行彻底的理论分析, 仅仅使用大量的随机模拟结果来说明模型或方法的有效性。当然, 因为模型都是有可变参数的, 随机模拟只能针对某些参数组合给出结果, 所以, 一般认为仅有模拟结果而没有理论分析结果的研究论文是不全面的。§3.5给出一个用随机模拟说明统计技术优良性的例子。

除了以上应用, 随机模拟还是许多新的统计方法的主要工具, 例如, 蒙特卡洛检验, bootstrap 置信区间和 bootstrap 偏差修正, MCMC。利用大量计算机计算 (包括随机模拟) 来进行统计推断的统计学分支叫做 “计算统计” (computational statistics), 在本章后面各节将介绍随机模拟的一些应用和技巧。

\section*{3.2 随机模拟积分}

某些非随机的问题也可以通过概率模型引入随机变量, 化为求随机模型的未知参数的问题。§3.1中用随机投点法估计 \(\pi\) 就是这样的一个例子。

随机模拟解决非随机问题的典型代表是计算定积分。通过对随机模拟定积分的讨论, 可以展示随机模拟中大部分的问题和技巧。随机模拟积分也称为蒙特卡洛 (Monte Carlo) 积分 (简称 MC 积分)。实际上, 统计中最常见的计算问题就是积分和最优化问题。

\section*{3.2.1 随机投点法}

设函数 \(h\left( x\right)\) 在有限区间 \(\left\lbrack  {a,b}\right\rbrack\) 上定义且有界,不妨设 \(0 \leq  h\left( x\right)  \leq  M\) 。要计算 \(I =\)  \({\int }_{a}^{b}h\left( x\right) {dx}\) ,相当于计算曲线下的区域 \(D = \{ \left( {x,y}\right)  : 0 \leq  y \leq  h\left( x\right) ,x \in  C = \left\lbrack  {a,b}\right\rbrack  \}\) 的面积。为此在 \(G = \left\lbrack  {a,b}\right\rbrack   \times  \left( {0,M}\right)\) 上均匀抽样 \(N\) 次,得随机点 \({Z}_{1},{Z}_{2}\ldots ,{Z}_{N},{Z}_{i} = \left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,N\) 。 令

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & {Z}_{i} \in  D \\  0, & \text{ 其它 } \end{array},\;i = 1,\ldots ,N}\right.
\]

则 \(\left\{  {\xi }_{i}\right\}\) 是独立重复试验结果, \(\left\{  {\xi }_{i}\right\}\) 独立同 \(\mathrm{b}\left( {1,p}\right)\) 分布,

\[
p = P\left( {{Z}_{i} \in  D}\right)  = V\left( D\right) /V\left( G\right)  = I/\left\lbrack  {M\left( {b - a}\right) }\right\rbrack   \tag{3.2}
\]

其中 \(V\left( \cdot \right)\) 表示区域面积。

\begin{center}
\includegraphics[max width=0.8\textwidth]{images/bo_d3dpgvk601uc738iknq0_3_311_866_1126_762_0.jpg}
\end{center}
\hspace*{3em} 

图 3.1: 随机投点法积分示意图

从模拟产生的随机样本 \({Z}_{1},{Z}_{2}\ldots ,{Z}_{N}\) ,可以用这 \(N\) 个点中落入曲线下方区域 \(D\) 的百分比 \(\widehat{p}\) 来估计 (3.2)中的概率 \(p\) (见图3.1),然后由 \(I = {pM}\left( {b - a}\right)\) 得到定积分 \(I\) 的近似值

\[
\widehat{I} = \widehat{p}M\left( {b - a}\right)  \tag{3.3}
\]

这种方法叫做随机投点法。这样计算的定积分有随机性, 误差中包含了随机模拟误差。

由强大数律可知

\[
\widehat{p} = \frac{\sum {\xi }_{i}}{N} \rightarrow  p,\text{ a.s. }\left( {N \rightarrow  \infty }\right)
\]

\[
\widehat{I} = \widehat{p}M\left( {b - a}\right)  \rightarrow  {pM}\left( {b - a}\right)  = I\text{ , a.s. }\left( {N \rightarrow  \infty }\right)
\]

即 \(N \rightarrow  \infty\) 时精度可以无限地提高 (当然,在计算机中要受到数值精度的限制)。

那么, 提高精度需要多大的代价呢? 由中心极限定理可知

\[
\sqrt{N}\left( {\widehat{p} - p}\right) /\sqrt{p\left( {1 - p}\right) }\overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,1}\right) ,\left( {N \rightarrow  \infty }\right) ,
\]

从而

\[
\sqrt{N}\left( {\widehat{I} - I}\right)  = M\left( {b - a}\right) \left( {\widehat{p} - p}\right) \overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,{\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}p\left( {1 - p}\right) }\right)  \tag{3.4}
\]

当 \(N\) 很大时 \(\widehat{I}\) 近似服从 \(\mathrm{N}\left( {I,{\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}p\left( {1 - p}\right) /N}\right)\) 分布,称此近似分布的方差 \(\lbrack M(b -\)  \(\left. {\left. a\right) {\rbrack }^{2}p\left( {1 - p}\right) /N}\right)\) 为 \(\widehat{I}\) 的渐近方差。计算渐近方差可以用 \(\widehat{p}\) 代替 \(p\) 估计为 \({\left\lbrack  M\left( b - a\right) \right\rbrack  }^{2}\widehat{p}\left( {1 - \widehat{p}}\right) /N\) 。 (3.4)说明 \(\widehat{I}\) 的误差为 \({O}_{p}\left( \frac{1}{\sqrt{N}}\right)\) ,这样,计算 \(\widehat{I}\) 的精度每增加一位小数,计算量需要增加 100 倍。随机模拟积分一般都服从这样的规律。

\section*{3.2.2 平均值法}

为了计算 \(I = {\int }_{a}^{b}h\left( x\right) {dx}\) ,上面用了类似于 \(§{2.2.4}\) 的舍选法的做法,在非随机问题中引入随机性时用了二维均匀分布和两点分布,靠求两点分布概率来估计积分 \(I\) 。随机投点法容易理解,但是效率较低,另一种效率更高的方法是利用期望值的估计。取 \(U \sim  \mathrm{U}\left( {a,b}\right)\) ,则

\[
E\left\lbrack  {h\left( U\right) }\right\rbrack   = {\int }_{a}^{b}h\left( u\right) \frac{1}{b - a}{du} = \frac{I}{b - a}
\]

\[
I = \left( {b - a}\right)  \cdot  {Eh}\left( U\right)
\]

若取 \(\left\{  {{U}_{i},i = 1,\ldots ,N}\right\}\) 独立同 \(\mathrm{U}\left( {a,b}\right)\) 分布,则 \({Y}_{i} = h\left( {U}_{i}\right) ,i = 1,2,\ldots ,N\) 是 iid 随机变量列, 由强大数律,

\[
\bar{Y} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right)  \rightarrow  {Eh}\left( U\right)  = \frac{I}{b - a},\;\text{ a.s. }\left( {N \rightarrow  \infty }\right)
\]

于是

\[
\widetilde{I} = \frac{b - a}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right)  \tag{3.5}
\]

是 \(I\) 的强相合估计。称这样计算定积分 \(I\) 的方法为平均值法。由中心极限定理有

\[
\sqrt{N}\left( {\widetilde{I} - I}\right) \overset{\mathrm{d}}{ \rightarrow  }\mathrm{N}\left( {0,{\left( b - a\right) }^{2}\operatorname{Var}\left( {h\left( U\right) }\right) }\right)
\]

其中

\[
\operatorname{Var}\left\lbrack  {h\left( U\right) }\right\rbrack   = {\int }_{a}^{b}{\left\lbrack  h\left( u\right)  - Eh\left( U\right) \right\rbrack  }^{2}\frac{1}{b - a}{du} \tag{3.6}
\]

仅与 \(h\) 有关,仍有 \(\widetilde{I} - I = {O}_{p}\left( \frac{1}{\sqrt{N}}\right)\) ,但是(3.5)的渐近方差小于(3.3)的渐近方差 (见 \(§{3.2.3}\) )。 \(\operatorname{Var}\left\lbrack  {h\left( U\right) }\right\rbrack\) 可以用模拟样本 \(\left\{  {{Y}_{i} = h\left( {U}_{i}\right) }\right\}\) 估计为

\[
\operatorname{Var}\left( {h\left( U\right) }\right)  \approx  \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {Y}_{i} - \bar{Y}\right) }^{2}. \tag{3.7}
\]

如果定积分区间是无穷区间,比如 \({\int }_{0}^{\infty }h\left( x\right) {dx}\) ,为了使用均匀分布随机数以及平均值法计算积分可以做积分变换,使积分区间变成有限区间。例如,作变换 \(t = 1/\left( {x + 1}\right)\) ,则

\[
{\int }_{0}^{\infty }h\left( x\right) {dx} = {\int }_{0}^{1}h\left( {\frac{1}{t} - 1}\right) \frac{1}{{t}^{2}}{dt}.
\]

从平均值法看出,定积分问题 \({\int }_{a}^{b}h\left( x\right) {dx}\) 等价于求 \({Eh}\left( U\right)\) ,其中 \(U \sim  \mathrm{U}\left( {a,b}\right)\) 。所以这一节讨论的方法也是用来求随机变量函数期望的随机模拟方法。对一般随机变量 \(X\) ,其取值范围不必局限于有限区间,为了求 \(X\) 的函数 \(h\left( X\right)\) 的期望 \(I = {Eh}\left( X\right)\) ,对 \(X\) 的随机数 \({X}_{i},i = 1,2,\ldots ,N\) ,令 \({Y}_{i} = h\left( {X}_{i}\right)\) ,也可以用平均值法 \(\widehat{I} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {X}_{i}\right)\) 来估计 \({Eh}\left( X\right) ,\widehat{I}\) 是 \(\operatorname{Eh}\left( X\right)\) 的无偏估计和强相合估计,若 \(\operatorname{Var}\left( {h\left( X\right) }\right)\) 存在,则 \(\widehat{I}\) 的方差为 \(\frac{1}{N}\operatorname{Var}\left( {h\left( X\right) }\right) ,\widehat{I}\) 有渐近正态分布 \(\mathrm{N}\left( {\operatorname{Eh}\left( X\right) ,\frac{1}{N}\operatorname{Var}\left( {h\left( X\right) }\right) }\right)\) 。设 \(\left\{  {Y}_{i}\right\}\) 的样本方差为 \({S}_{N}^{2}\) ,可以用 \({S}_{N}^{2}/N\) 估计 \(\widehat{I}\) 的方差,用 \(\widehat{I} \pm  2{S}_{N}/\sqrt{N}\) 作为 \(I\) 的近似 \({95}\%\) 置信区间。

例 3.2.1. 设 \(X \sim  \mathrm{N}\left( {0,1}\right)\) ,求 \(I = E{\left| X\right| }^{\frac{3}{2}}\) 。

解: 作变量替换积分可得

\[
I = 2{\int }_{0}^{\infty }{x}^{\frac{3}{2}}\frac{1}{\sqrt{2\pi }}{e}^{-\frac{1}{2}{x}^{2}}{dx}\;\left( {\text{ 令 }t = \frac{1}{2}{x}^{2}}\right)
\]

\[
= \frac{{2}^{\frac{3}{4}}\Gamma \left( \frac{5}{4}\right) }{\sqrt{\pi }} \approx  {0.86004}
\]

如果用平均值法估计 \(I\) ,取抽样样本量 \(N = {10000}\) ,产生标准正态随机数 \({X}_{i},i =\)  \(1,2,\ldots ,n\) ,令 \({Y}_{i} = {\left| {X}_{i}\right| }^{\frac{3}{2}}\) ,令 \(\widehat{I} = \bar{Y} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left| {X}_{i}\right| }^{\frac{3}{2}},{S}_{N}^{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\left( {\left| {X}_{i}\right| }^{\frac{3}{2}} - \widehat{I}\right) }^{2}\) ,一次模拟的结果得到 \(\widehat{I} = {0.8658},{S}_{N} = {0.9409}\) ,渐近标准差为 \({S}_{N}/\sqrt{N} = {0.00941}\) ,说明 \({95}\%\) 置信水平下误差界为 \(\pm  {0.02}\) ,结果只有将近两位小数的精度。为了达到四位小数的精度,需要误差控制在 \(\pm  {0.00005}\) 以下,需要 \(2{S}_{N}/\sqrt{N} < {0.00005}\) ,代入得 \(N\) 需要超过 90 亿,可见随机模拟方法提高精度的困难程度。

以上的评估误差的方法使用了渐近方差和渐近正态分布, 在有些模拟问题中估计量不一定有渐近方差, 或渐近分布不一定准确 (比如产生的随机数不独立时)。在计算量允许的情形下,可以用不同的随机数种子重复得到 \(B\) 个 \(\widehat{I}\) 的估计值,用这 \(B\) 个 \(\widehat{I}\) 的样本标准差来估计 \(\widehat{I}\) 的抽样分布标准差。取 \(B = {100}\) 的一次计算结果得到的 \(\widehat{I}\) 的抽样分布标准差为 0.00924, 与从一次模拟得到的渐近标准差结果十分接近。也可以使用 \(§{3.6}\) 的 bootstrap 方法获取 \(B\) 组 bootstrap 样本,得到 \(B\) 个 bootstrap 的 \(\widehat{I}\) 样本值,从中估计抽样分布标准差。bootstrap 方法可以避免重新计算 \(h\left( {X}_{i}\right)\) 的值,在更复杂的问题中往往 \(h\left( {X}_{i}\right)\) 的计算是速度很慢的,比如说, \(h\left( {X}_{i}\right)\) 本身也需要用随机模拟或者数值优化方法计算。

实际上,一维积分用数值方法均匀布点计算一般更有效。比如,令 \({x}_{i} = a + \frac{b - a}{N}i\) , \(i = 0,1,\ldots ,N\) ,估计 \(I\) 为 (复合梯形求积公式)

\[
{I}_{0} = \frac{b - a}{N}\left\lbrack  {\frac{1}{2}h\left( a\right)  + \mathop{\sum }\limits_{{i = 1}}^{{N - 1}}h\left( {x}_{i}\right)  + \frac{1}{2}h\left( b\right) }\right\rbrack   \tag{3.8}
\]

则当 \(h\left( x\right)\) 在 \(\left\lbrack  {a,b}\right\rbrack\) 上二阶连续可微时误差 \({I}_{0} - I\) 仅为 \(O\left( \frac{1}{{N}^{2}}\right)\) 阶,比随机模拟方法得到的精度要高得多,而且当 \(h\left( x\right)\) 有更好的光滑性时还可以用更精确的求积公式得到更高精度。所以, 被积函数比较光滑的一元定积分问题一般不需要用随机模拟来计算。

\section*{3.2.3 高维定积分}

上面的两种计算一元函数定积分的方法可以很容易地推广到多元函数定积分, 或称高维定积分。设 \(d\) 元函数 \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)\) 定义于超矩形

\[
C = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  : {a}_{i} \leq  {x}_{i} \leq  {b}_{i},i = 1,2,\ldots ,d}\right\}
\]

且

\[
0 \leq  h\left( {{x}_{1},\ldots ,{x}_{d}}\right)  \leq  M,\forall x \in  C.
\]

令

\[
D = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d},y}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  \in  C,0 \leq  y \leq  h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right) }\right\}  ,
\]

\[
G = \left\{  {\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d},y}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  \in  C,0 \leq  y \leq  M}\right\}
\]

为计算 \(d\) 维定积分

\[
I = {\int }_{{a}_{d}}^{{b}_{d}}\cdots {\int }_{{a}_{2}}^{{b}_{2}}{\int }_{{a}_{1}}^{{b}_{1}}h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right) d{x}_{1}d{x}_{2}\cdots d{x}_{d}, \tag{3.9}
\]

产生服从 \(d + 1\) 维空间中的超矩形 \(G\) 内的均匀分布的独立抽样 \({\mathbf{Z}}_{1},{\mathbf{Z}}_{2},\ldots ,{\mathbf{Z}}_{N}\) ,令

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & {\mathbf{Z}}_{i} \in  D \\  0, & {\mathbf{Z}}_{i} \in  G - D \end{array},\;i = 1,2,\ldots ,N}\right.
\]

则 \({\xi }_{i}\) iid \(\mathrm{b}\left( {1,p}\right)\) ,

\[
p = P\left( {{\mathbf{Z}}_{i} \in  D}\right)  = \frac{V\left( D\right) }{V\left( G\right) } = \frac{I}{{MV}\left( C\right) } = \frac{I}{M\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) }
\]

其中 \(V\left( \cdot \right)\) 表示区域体积。令 \(\widehat{p}\) 为 \(N\) 个随机点中落入 \(D\) 的百分比,则

\[
\widehat{p} = \frac{\sum {\xi }_{i}}{N} \rightarrow  p,\text{ a.s. }\left( {N \rightarrow  \infty }\right) ,
\]

用

\[
{\widehat{I}}_{1} = \widehat{p}V\left( G\right)  = \widehat{p} \cdot  {MV}\left( C\right)  = \widehat{p} \cdot  M\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right)  \tag{3.10}
\]

估计积分 \(I\) ,则 \({\widehat{I}}_{1}\) 是 \(I\) 的无偏估计和强相合估计。称用(3.10)计算高维定积分 \(I\) 的方法为随机投点法。由中心极限定理知

\[
\sqrt{N}\left( {\widehat{p} - p}\right) /\sqrt{p\left( {1 - p}\right) }\overset{d}{ \rightarrow  }\mathrm{\;N}\left( {0,1}\right)
\]

\[
\sqrt{N}\left( {{\widehat{I}}_{1} - I}\right) \overset{d}{ \rightarrow  }\mathrm{\;N}\left( {0,{\left( M\mathop{\prod }\limits_{{j = 1}}^{n}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}p\left( {1 - p}\right) }\right) ,
\]

\({\widehat{I}}_{1}\) 的渐近方差为

\[
\frac{{\left( M\mathop{\prod }\limits_{{j = 1}}^{n}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}p\left( {1 - p}\right) }{N} \tag{3.11}
\]

所以 \({\widehat{I}}_{1}\) 的随机误差仍为 \({O}_{p}\left( \frac{1}{\sqrt{N}}\right) ,N \rightarrow  \infty\) 时的误差阶不受维数 \(d\) 的影响,这是随机模拟方法与其它数值计算方法相比一个重大优势。

在计算高维积分(3.9)时,仍可以通过估计 \({Eh}\left( \mathbf{U}\right)\) 来获得,其中 \(\mathbf{U}\) 服从 \({R}^{d}\) 中的超矩形 \(C\) 上的均匀分布。设 \({\mathbf{U}}_{i} \sim\) iid \(\mathrm{U}\left( C\right) ,i = 1,2,\ldots ,N\) ,则 \(h\left( {\mathbf{U}}_{i}\right) ,i = 1,2,\ldots ,N\) 是 iid 随机变量列,

\[
{Eh}\left( {\mathbf{U}}_{i}\right)  = {\int }_{C}h\left( \mathbf{u}\right) \frac{1}{\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) }d\mathbf{u} = \frac{I}{\mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) },
\]

估计 \(I\) 为

\[
{\widehat{I}}_{2} = \mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right)  \cdot  \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {U}_{i}\right) , \tag{3.12}
\]

称用(3.12)计算高维定积分 \(I\) 的方法为平均值法。由强大数律

\[
{\widehat{I}}_{2} \rightarrow  \mathop{\prod }\limits_{{j = 1}}^{d}\left( {{b}_{j} - {a}_{j}}\right) {Eh}\left( U\right)  = I,\text{ a.s. }\left( {N \rightarrow  \infty }\right) ,
\]

\({\widehat{I}}_{2}\) 的渐近方差为

\[
\frac{{\left( V\left( C\right) \right) }^{2}\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right) }{N} = \frac{{\left( \mathop{\prod }\limits_{{j = 1}}^{d}\left( {b}_{j} - {a}_{j}\right) \right) }^{2}\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right) }{N}. \tag{3.13}
\]

\(N \rightarrow  \infty\) 时的误差阶也不受维数 \(d\) 的影响。

我们来比较随机投点法(3.10)与平均值法(3.12)的精度,只要比较其渐近方差。对 \(I =\)  \({\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) ,设 \({\widehat{I}}_{1}\) 为随机投点法的估计, \({\widehat{I}}_{2}\) 为平均值法的估计。因设 \(0 \leq  h\left( \mathbf{x}\right)  \leq  M\) ,不妨设 \(0 \leq  h\left( \mathbf{x}\right)  \leq  1\) ,取 \(h\left( \mathbf{x}\right)\) 的上界 \(M = 1\) 。

令 \({\mathbf{X}}_{i} \sim\) iid \(\mathrm{U}\left( C\right) ,{\eta }_{i} = h\left( {\mathbf{X}}_{i}\right) ,{Y}_{i} \sim\) iid \(\mathrm{U}\left( {0,1}\right)\) 与 \(\left\{  {\mathbf{X}}_{i}\right\}\) 独立,

\[
{\xi }_{i} = \left\{  {\begin{array}{ll} 1, & \text{ 当 }{Y}_{i} \leq  h\left( {\mathbf{X}}_{i}\right) \\  0, & \text{ 当 }{Y}_{i} > h\left( {\mathbf{X}}_{i}\right)  \end{array}i = 1,2,\ldots ,N,}\right.
\]

这时有

\[
{\widehat{I}}_{1} = V\left( C\right) \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\xi }_{i},\;{\widehat{I}}_{2} = V\left( C\right) \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\eta }_{i}
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{1}{N}{V}^{2}\left( C\right)  \cdot  \frac{I}{V\left( C\right) }\left( {1 - \frac{I}{V\left( C\right) }}\right)
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{1}{N}{V}^{2}\left( C\right)  \cdot  \left( {\frac{1}{V\left( C\right) }{\int }_{C}{h}^{2}\left( \mathbf{x}\right) d\mathbf{x} - {\left( \frac{I}{V\left( C\right) }\right) }^{2}}\right)
\]

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  - \operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{V\left( C\right) }{N}{\int }_{C}\left\{  {h\left( \mathbf{x}\right)  - {h}^{2}\left( \mathbf{x}\right) }\right\}  d\mathbf{x} \geq  0
\]

可见平均值法精度更高。事实上,随机投点法多用了随机数 \({Y}_{i}\) ,必然会增加抽样随机误差。

在计算高维积分(3.10)时,如果用网格法作数值积分,把超矩形 \(C = \left\lbrack  {{a}_{1},{b}_{1}}\right\rbrack   \times  \left\lbrack  {{a}_{2},{b}_{2}}\right\rbrack   \times\)  \(\cdots  \times  \left\lbrack  {{a}_{d},{b}_{d}}\right\rbrack\) 的每个边分成 \(n\) 个格子,就需要 \(N = {n}^{d}\) 个格子点,如果用每个小超矩形的中心作为代表点,可以达到 \(O\left( {n}^{-2}\right)\) 精度,即 \(O\left( {N}^{-2/d}\right)\) ,当维数增加时为了提高一倍精度需要 \({2}^{d/2}\) 倍的代表点。比如 \(d = 8\) ,精度只有 \(O\left( {N}^{-1/4}\right)\) 。高维的问题当维数增加时计算量会迅猛增加,以至于无法计算,这个问题称为维数诅咒。如果用 Monte Carlo 积分,则精度为 \({O}_{p}\left( {N}^{-1/2}\right)\) ,与 \(d\) 关系不大。所以 Monte Carlo 方法在高维积分中有重要应用。为了提高积分计算精度,需要减小 \({O}_{\mathrm{p}}\left( {N}^{-1/2}\right)\) 中的常数项,即减小 \(\widehat{I}\) 的渐近方差。

例 3.2.2. 考虑

\[
f\left( {{x}_{1},{x}_{2},\ldots ,{x}_{d}}\right)  = {x}_{1}^{2}{x}_{2}^{2}\cdots {x}_{d}^{2},0 \leq  {x}_{1} \leq  1,0 \leq  {x}_{2} \leq  1,\ldots ,0 \leq  {x}_{d} \leq  1
\]

的积分

\[
I = {\int }_{0}^{1}\cdots {\int }_{0}^{1}{\int }_{0}^{1}{x}_{1}^{2}{x}_{2}^{2}\cdots {x}_{d}^{2}d{x}_{1}d{x}_{2}\cdots d{x}_{d}
\]

当然,这个积分是有精确解 \(I = {\left( 1/3\right) }^{d}\) 的,对估计 \(I\) 的结果我们可以直接计算误差。以 \(d = 8\) 为例比较以下三种方法的精度,这时真值 \(I = {\left( 1/3\right) }^{8} \approx  {1.524} \times  {10}^{-4}\) 。

用 \(n\) 网格点法, \(N = {n}^{d}\) ,公式为

\[
{I}_{0} = \frac{1}{N}\mathop{\sum }\limits_{{{i}_{1} = 1}}^{n}\mathop{\sum }\limits_{{{i}_{2} = 1}}^{n}\ldots \mathop{\sum }\limits_{{{i}_{d} = 1}}^{n}f\left( {\frac{2{i}_{1} - 1}{2n},\frac{2{i}_{2} - 1}{2n},\ldots ,\frac{2{i}_{d} - 1}{2n}}\right)
\]

误差绝对值为 \({e}_{0} = \left| {{I}_{0} - I}\right|\) 。如果取 \(n = 5,d = 8\) ,需要计算 \(N = {390625}\) 个点的函数值,计算量相当大。

用随机投点法求 \(I\) ,先在 \({\left( 0,1\right) }^{d} \times  \left( {0,1}\right)\) 均匀抽样 \(\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\ldots ,{\xi }_{i}^{\left( d\right) },{y}_{i}}\right) ,i = 1,\ldots ,N\) ,令 \({\widehat{I}}_{1}\) 为 \({y}_{i} \leq  f\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\ldots ,{\xi }_{i}^{\left( d\right) }}\right)\) 成立的百分比。因为 \({\widehat{I}}_{1}\) 是随机的,误差绝对值 \(\left| {{\widehat{I}}_{1} - I}\right|\) 也是随机的,所以我们重复试验 \(B\) 次,计算 \(B\) 次的误差绝对值的平均值 \({e}_{1}\) ,作为 \(E\left| {{\widehat{I}}_{1} - I}\right|\) 的估计值。取 \(B\) 多大合适呢? 因为计算量很大,先取了 \({B}_{1} = {10}\) ,用得到的 \({B}_{1}\) 个 \(\left| {{\widehat{I}}_{1} - I}\right|\) 样本标准差 \({s}_{1}\) 可以估计 \(B\) 个 \(\left| {{\widehat{I}}_{1} - I}\right|\) 的样本平均值的标准差为 \({\mathrm{{SE}}}_{1} = {s}_{1}/\sqrt{B}\) ,发现 \(B = {1000}\) 时可以控制在相对误差 2\% 以下。

用平均值方法, 公式为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}f\left( {{\xi }_{i}^{\left( 1\right) },{\xi }_{i}^{\left( 2\right) },\cdots ,{\xi }_{i}^{\left( d\right) }}\right)
\]

其中 \({\xi }_{i}^{\left( j\right) },i = 1,\ldots ,N,j = 1,\ldots ,d\) 独立同 \(\mathrm{U}\left( {0,1}\right)\) 分布。类似地,重复 \(B = {1000}\) 次,计算 \(B\) 次的误差绝对值 \(\left| {{\widehat{I}}_{2} - I}\right|\) 的平均值 \({e}_{2}\) ,作为 \(E\left| {{\widehat{I}}_{2} - I}\right|\) 的估计。

最后,三种方法计算量相同(都计算 \(N = {5}^{8}\) 次函数值)的情况下,得到网格点法的误差绝对值为 \({e}_{0} = {1.2} \times  {10}^{-5}\) ,随机投点法的误差绝对值平均值为 \({e}_{1} = {1.6} \times  {10}^{-5}\) ,平均值法的误差绝对值平均值为 \({e}_{2} = {0.20} \times  {10}^{-5}\) ,此问题结果中平均值法比其它两种方法精度高了至少 5 倍。

\section*{3.2.4 重要抽样法}

\(I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) 中积分区域 \(C\) 可能是任意形状的,也可能无界, \(h\left( \mathbf{x}\right)\) 在 \(C\) 内各处的取值大小差异可能很大,使得直接用平均值法估计 \(I\) 时,很多样本点处于 \(\left| {h\left( \mathbf{x}\right) }\right|\) 接近于零的地方, 造成浪费,另外使得 \({\widehat{I}}_{2}\) 的渐近方差 (见(3.13)) 中的 \(\operatorname{Var}\left( {h\left( \mathbf{U}\right) }\right)\) 很大 \(\left( {\mathbf{U} \sim  \mathrm{U}\left( C\right) }\right)\) 。为此,考虑用非均匀抽样: \(\left| {h\left( \mathbf{x}\right) }\right|\) 大的地方密集投点, \(\left| {h\left( \mathbf{x}\right) }\right|\) 小的地方稀疏投点。这样可以有效利用样本。设 \(g\left( \mathbf{x}\right) ,\mathbf{x} \in  C\) 是一个密度,其形状与 \(\left| {h\left( \mathbf{x}\right) }\right|\) 相近,且当 \(g\left( \mathbf{x}\right)  = 0\) 时 \(h\left( \mathbf{x}\right)  = 0\) ,当 \(\parallel \mathbf{x}\parallel  \rightarrow  \infty\) 时 \(h\left( \mathbf{x}\right)  = o\left( {g\left( \mathbf{x}\right) }\right)\) 。称 \(g\left( \mathbf{x}\right)\) 为试投密度或重要抽样密度。

设 \({\mathbf{X}}_{i}\) iid \(\sim  g\left( \mathbf{x}\right) ,i = 1,2,\ldots ,N\) 。令

\[
{\eta }_{i} = \frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) },i = 1,2,\ldots ,N
\]

则

\[
E{\eta }_{1} = {\int }_{C}\frac{h\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) }g\left( \mathbf{x}\right) d\mathbf{x} = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x} = I
\]

因此可以用 \(\left\{  {{\eta }_{i},i = 1,2,\ldots ,N}\right\}\) 的样本平均值来估计 \(I\) ,即

\[
{\widehat{I}}_{3} = \bar{\eta } = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{h\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }. \tag{3.14}
\]

\({\widehat{I}}_{3}\) 是 \(I\) 的无偏估计和强相合估计。 \({\widehat{I}}_{3}\) 的渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right) \frac{1}{N} = O\left( \frac{1}{N}\right) , \tag{3.15}
\]

当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\frac{\left| h\left( \mathbf{x}\right) \right| }{g\left( \mathbf{x}\right) }\) 近似为常数,方差 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小,这时 \({\widehat{I}}_{3}\) 的随机误差可以很小。用(3.14)估计 \(I\) 与 \(§{2.2.4}\) 的舍选法 II 有类似的想法,这种方法叫做重要抽样法 (importance sampling), 是随机模拟的重要方法。

为什么当 \(g\left( \mathbf{x}\right)\) 和 \(\left| {h\left( \mathbf{x}\right) }\right|\) 形状很接近时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 很小？事实上,

\[
\operatorname{Var}\left( \frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - \left( {E\frac{h\left( \mathbf{X}\right) }{g\left( \mathbf{X}\right) }}\right)  = E\left( \frac{{h}^{2}\left( \mathbf{X}\right) }{{g}^{2}\left( \mathbf{X}\right) }\right)  - {\left( {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\right) }^{2}
\]

\[
\geq  {\left( E\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\right) }^{2} - {I}^{2}\;\text{ (由 Jensen 不等式) } \tag{3.16}
\]

\[
= {\left( {\int }_{C}\left| h\left( \mathbf{x}\right) \right| d\mathbf{x}\right) }^{2} - {I}^{2}, \tag{3.17}
\]

当且仅当 \(\frac{\left| h\left( \mathbf{X}\right) \right| }{g\left( \mathbf{X}\right) }\) 为常数时不等式(3.16)中的等号成立,即 \(g\left( \mathbf{x}\right)  = \frac{1}{{\int }_{C}\left| {h\left( \mathbf{t}\right) }\right| {dt}}\left| {h\left( \mathbf{x}\right) }\right| ,\mathbf{x} \in  C\) 时 \(\operatorname{Var}\left( {\widehat{I}}_{3}\right)\) 达到最小。

上述求积分的问题也可以表述为求随机变量函数的期望问题。设 \(\mathbf{Y} \sim  f\left( \mathbf{y}\right)\) ,要求 \(\mathbf{Y}\) 的函数 \(h\left( \mathbf{Y}\right)\) 的期望值 \({Eh}\left( \mathbf{Y}\right)  = \int h\left( \mathbf{y}\right) f\left( \mathbf{y}\right) d\mathbf{y}\) ,可以抽取 \(\mathbf{Y}\) 的随机数 \({\mathbf{Y}}_{i},i = 1,2,\ldots ,N\) ,然后用平均值法 \(\frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{Y}}_{i}\right)\) 估计 \({Eh}\left( \mathbf{Y}\right)\) 。但是,如果很难生成 \(\mathbf{Y}\) 的随机数,或者 \(\mathbf{Y}\) 的分布集中于 \(h\left( \mathbf{y}\right)\) 接近于零的位置以至于积分效率很低,可以找一个试投密度 \(g\left( \mathbf{x}\right)\) ,从 \(g\left( \mathbf{x}\right)\) 产生随机数 \({\mathbf{X}}_{i},i = 1,2,\ldots ,N\) ,用如下重要抽样法

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}h\left( {\mathbf{X}}_{i}\right) \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) } \tag{3.18}
\]

估计 \({Eh}\left( \mathbf{Y}\right)\) ,易见 \(E{\widehat{I}}_{3.1} = {Eh}\left( \mathbf{Y}\right)\) 。称 \({W}_{i} = \frac{f\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 为重要性权重, \({Eh}\left( \mathbf{Y}\right)\) 可以用重要性权重估计为

\[
{\widehat{I}}_{3.1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) . \tag{3.19}
\]

选取试投密度 \(g\left( \mathbf{x}\right)\) 时,要求当 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  \neq  0\) 时 \(g\left( \mathbf{x}\right)  \neq  0\) ,当 \(\mathbf{x}\) 趋于无穷时 \(h\left( \mathbf{x}\right) f\left( \mathbf{x}\right)  =\)  \(o\left( {g\left( \mathbf{x}\right) }\right)\) ( \(g\left( \mathbf{x}\right)\) 相对厚尾),一般还要求 \(g\left( \mathbf{x}\right)\) 的形状与 \(\left| {h\left( \mathbf{x}\right) }\right| f\left( \mathbf{x}\right)\) 形状接近。如果 \(g\left( \mathbf{x}\right)\) 的相对厚尾性难以确定, 可以使用如下保险的试投密度

\[
\widetilde{g}\left( \mathbf{x}\right)  = {\rho g}\left( \mathbf{x}\right)  + \left( {1 - \rho }\right) r\left( \mathbf{x}\right) , \tag{3.20}
\]

其中 \(\rho\) 接近于 1, \(r\left( \mathbf{x}\right)\) 是柯西分布或 Pareto 分布这样的重尾分布。要生成 \(N\) 个 \(\widetilde{g}\left( \mathbf{x}\right)\) 的随机数,只要生成 \({N\rho }\) 个 \(g\left( \mathbf{x}\right)\) 的随机数和 \(N\left( {1 - \rho }\right)\) 个 \(r\left( \mathbf{x}\right)\) 的随机数。

选取不适当的试投密度会把绝大多数样本点投到了对计算积分不重要的位置, 使得样本点中只有极少数点是真正有作用的。在多维问题中合适的试投密度尤其难找, 经常需要反复试验。

有时 \(\mathbf{Y} \sim  f\left( \cdot \right)\) 不仅很难直接抽样,而且 \(f\left( \cdot \right)\) 本身未知,只能确定到差一个常数倍的 \(\widetilde{f}\left( \mathbf{x}\right)  = {cf}\left( \mathbf{x}\right)\) ,常数 \(c\) 未知,为了求常数 \(c\) 需要计算 \(c = \int \widetilde{f}\left( \mathbf{y}\right) d\mathbf{y}\) ,计算 \(c\) 一般很困难。这时,定义重要性权重为 \({W}_{i} = \frac{\widetilde{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) ,公式(3.19)可以改成

\[
{\widehat{I}}_{4} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}} \tag{3.21}
\]

这称为标准化重要抽样法。对(3.21)的分子和分母都除以 \({cN}\) 后分子 a.s. 收敛到 \({Eh}\left( \mathbf{Y}\right)\) ,分母 a.s. 收敛到 1,所以 (3.21) 是 \(\operatorname{Eh}\left( \mathbf{Y}\right)\) 的强相合估计,但不是无偏的。标准化重要抽样估计往往比无偏估计 \({\widehat{I}}_{3.1}\) 有更小的均方误差。关于标准化重要抽样法的渐近方差的讨论参见 \(\operatorname{Liu}{\left( {2001}\right) }^{\left\lbrack  {28}\right\rbrack  }§{2.5.3}\) 。

如果需要对多个不同的函数 \(h\left( \cdot \right)\) 计算 \({Eh}\left( \mathbf{Y}\right)\) ,则选取试抽样密度 \(g\left( \mathbf{x}\right)\) 时应使得 \(g\left( \mathbf{x}\right)\) 尽可能与 \(\mathbf{Y}\) 的密度 \(f\left( \mathbf{x}\right)\) 形状接近,这样权重 \({W}_{i} = \frac{\widehat{f}\left( {\mathbf{X}}_{i}\right) }{g\left( {\mathbf{X}}_{i}\right) }\) 的分布不过于偏斜,不至于出现绝大部分权重集中于少数样本点的情形。抽样值 \({\mathbf{X}}_{i}\) 与权重 \({W}_{i}\) 一起可以看作是分布 \(f\left( \cdot \right)\) 的某种抽样。

定义 (适当加权抽样) 随机变量序列 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}}\right) ,i = 1,2,\ldots ,N}\right\}\) 称为关于密度 \(f\left( \cdot \right)\) 的适当加权抽样,如果对于任何平方可积函数 \(h\left( \cdot \right)\) 都有

\[
E\left\lbrack  {h\left( {X}_{i}\right) {W}_{i}}\right\rbrack   = c{E}_{f}\left\lbrack  {h\left( X\right) }\right\rbrack   = c\int h\left( \mathbf{x}\right) f\left( \mathbf{x}\right) d\mathbf{x},i = 1,2,\ldots ,N,
\]

其中 \(c\) 是归一化常数。

设随机变量(X, W)联合密度为 \(g\left( {\mathbf{x},w}\right)\) ,则(X, W)的样本为密度 \(f\left( \cdot \right)\) 的适当加权抽样的充分必要条件是

\[
{E}_{g}\left( {W \mid  \mathbf{x}}\right)  = {E}_{g}\left( W\right) \frac{f\left( \mathbf{x}\right) }{g\left( \mathbf{x}\right) },\forall \mathbf{x},
\]

其中 \({E}_{g}\left( W\right)\) 是关于 \(W\) 的边缘密度的期望, \({E}_{g}\left( {W \mid  \mathbf{x}}\right)\) 是在(X, W)的联合密度下条件期望 \(E\left( {W \mid  \mathbf{X}}\right)\) 在 \(\mathbf{X} = \mathbf{x}\) 处的值。

在重要抽样法和标准化重要抽样法的实际应用中, 好的试抽样分布很难获得, 所以权重 \(\left\{  {{W}_{i} = f\left( {\mathbf{X}}_{i}\right) /g\left( {\mathbf{X}}_{i}\right) }\right\}\) 经常会差别很大,使得抽样样本主要集中在少数几个权重最大的样本点上。为此, 可以舍弃权重太小的样本点, 重新抽样替换这样的样本点。这种方法称为舍选控制 (Rejection Control),描述如下。

首先,需要选定权重的一个阈值 \(c\) ,然后对每个样本点 \({\mathbf{X}}_{i}\) ,若 \({W}_{i} \geq  c\) ,则接受 \({\mathbf{X}}_{i}\) ; 若 \({W}_{i} < c\) ,则以概率 \({W}_{i}/c\) 接受 \({\mathbf{X}}_{i}\) ,否则舍弃 \({\mathbf{X}}_{i}\) ,重新抽取。最后,所有样本点的权重调整为新的

\[
{W}_{i}^{ * } = {p}_{c}\frac{{W}_{i}}{\min \left\{  {1,{W}_{i}/c}\right\}  } = {p}_{c}\max \left\{  {{W}_{i},c}\right\}  , \tag{3.22}
\]

其中

\[
{p}_{c} = \int \min \left\{  {1,\frac{f\left( \mathbf{x}\right) }{{cg}\left( \mathbf{x}\right) }}\right\}  g\left( \mathbf{x}\right) d\mathbf{x} = \int \min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  d\mathbf{x} \tag{3.23}
\]

是归一化常数,如果使用标准化重要抽样法, \({p}_{c}\) 可以省略,否则, \({p}_{c}\) 可以估计为

\[
{\widehat{p}}_{c} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\min \left\{  {1,\frac{{W}_{i}}{c}}\right\}  . \tag{3.24}
\]

这样得到的 \(\left\{  {\left( {{\mathbf{X}}_{i},{W}_{i}^{ * }}\right) ,i = 1,\ldots ,N}\right\}\) 是关于 \(f\left( \cdot \right)\) 适当加权的,被接受的 \({\mathbf{X}}_{i}\) 的分布密度为

\[
{g}^{ * }\left( \mathbf{x}\right)  = \frac{1}{{p}_{c}}\min \left\{  {g\left( \mathbf{x}\right) ,\frac{f\left( \mathbf{x}\right) }{c}}\right\}  . \tag{3.25}
\]

阈值 \(c\) 在实际中可以从权重 \(\left\{  {W}_{i}\right\}\) 选取,比如,取为 \(\left\{  {W}_{i}\right\}\) 的某个分位数。

例 3.2.3. 用 \(\mathrm{{MC}}\) 积分法计算 \(I = {\int }_{0}^{1}{e}^{x}{dx} = e - 1 \approx  {1.718}\) 。对被积函数 \(h\left( x\right)  = {e}^{x}\) 做泰勒展开得

\[
{e}^{x} = 1 + x + \frac{{x}^{2}}{2!} + \cdots
\]

取

\[
g\left( x\right)  = c\left( {1 + x}\right)  = \frac{2}{3}\left( {1 + x}\right)
\]

要产生 \(g\left( x\right)\) 的随机数可以用逆变换法,密度 \(g\left( x\right)\) 的分布函数 \(G\left( x\right)\) 的反函数为

\[
{G}^{-1}\left( y\right)  = \sqrt{1 + {3y}} - 1,0 < y < 1
\]

因此,取 \({U}_{i}\) iid \(\mathrm{U}\left( {0,1}\right)\) ,令 \({X}_{i} = \sqrt{1 + 3{U}_{i}} - 1,i = 1,2,\ldots ,N\) ,则重要抽样法的积分公式为

\[
{\widehat{I}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{{e}^{{X}_{i}}}{\frac{2}{3}\left( {1 + {X}_{i}}\right) }
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \frac{1}{N}\left( {\frac{3}{2}{\int }_{0}^{1}\frac{{e}^{2x}}{1 + x}{dx} - {I}^{2}}\right)  \approx  {0.02691}/N.
\]

如果用平均值法, 估计公式为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{e}^{{U}_{i}}
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{1}{N}\left( {{\int }_{0}^{1}{e}^{2x}{dx} - {I}^{2}}\right)  \approx  {0.2420}/N \approx  {9.0} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

是重要抽样法方差的 9 倍。

如果用随机投点法, \(h\left( x\right)  = {e}^{x} \leq  e\left( {0 < x < 1}\right)\) ,取上界 \(M = e\) ,向 \(\left\lbrack  {0,1}\right\rbrack   \times  \left\lbrack  {0,M}\right\rbrack\) 随机投点,落到 \(f\left( x\right)\) 下方的概率为

\[
p = I/\left( {M\left( {b - a}\right) }\right)  = \left( {e - 1}\right) /e,
\]

设投 \(N\) 点落到 \(h\left( x\right)\) 下方的频率为 \(\widehat{p}\) ,用随机投点法估计 \(I\) 的公式为

\[
{\widehat{I}}_{1} = \widehat{p} \cdot  M\left( {b - a}\right)  = e\widehat{p},
\]

渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = {e}^{2}p\left( {1 - p}\right) /N = \left( {e - 1}\right) /N \approx  {1.718}/N \approx  {7.1} \times  \operatorname{Var}\left( {\widehat{I}}_{2}\right)  \approx  {64.8} \times  \operatorname{Var}\left( {\widehat{I}}_{3}\right)
\]

可见选择合适的抽样算法对减少计算量、提高精度是十分重要的。

例 3.2.4. 设二元函数 \(f\left( {x,y}\right)\) 定义如下

\[
f\left( {x,y}\right)  = \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {45}{\left( y + {0.1}\right) }^{4}}\right\}
\]

求如下二重定积分

\[
I = {\int }_{-1}^{1}{\int }_{-1}^{1}f\left( {x,y}\right) {dxdy}
\]

\(f\left( {x,y}\right)\) 有两个分别以(-0.4,0.5)和(0.5, - 0.1)为中心的峰,对积分有贡献的区域主要集中在 (-0.4,0.5)和(0.5, - 0.1)附近,在其他地方函数值很小,对积分贡献很小。

用平均值法 (3.12),取点数 \(N = {10000}\) 时 \({\widehat{I}}_{2}\) 的一个估计值为 \({\widehat{I}}_{2} = {0.132}\) ,从这一次模拟估计的 \({\widehat{I}}_{2}\) 的渐近方差为 \({1.86} \times  {10}^{-5}\) 。重复模拟 \(B = {100}\) 次,得到这 100 个 \({\widehat{I}}_{2}\) 的平均值为 0.1256,样本方差为 \({1.95} \times  {10}^{-5}\) 。

用重要抽样法。取试投密度为

\[
g\left( {x,y}\right)  \propto  \widetilde{g}\left( {x,y}\right)
\]

\[
= \exp \left\{  {-{45}{\left( x + {0.4}\right) }^{2} - {60}{\left( y - {0.5}\right) }^{2}}\right\}   + {0.5}\exp \left\{  {-{90}{\left( x - {0.5}\right) }^{2} - {10}{\left( y + {0.1}\right) }^{2}}\right\}  ,
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

这样抽取到 \(\left\lbrack  {-1,1}\right\rbrack   \times  \left\lbrack  {-1,1}\right\rbrack\) 范围外的点对积分没有贡献,因为构成 \(g\left( {x,y}\right)\) 的两个密度都很集中,所以效率损失不大。需要求使得 \(\widetilde{g}\left( {x,y}\right)\) 化为密度的比例常数。记 \(\mathrm{N}\left( {\mu ,{\sigma }^{2}}\right)\) 的分布密度为 \(f\left( {x;\mu ,{\sigma }^{2}}\right)\) ,对 \(\widetilde{g}\left( {x,y}\right)\) 积分得

\[
{\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\widetilde{g}\left( {x,y}\right)  = \sqrt{{2\pi }/{90}}{\int }_{-\infty }^{\infty }f\left( {x; - {0.4},{90}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{120}}{\int }_{-\infty }^{\infty }f\left( {y;{0.5},{120}^{-1}}\right) {dy}
\]

\[
+ {0.5}\sqrt{{2\pi }/{180}}{\int }_{-\infty }^{\infty }f\left( {x;{0.5},{180}^{-1}}\right) {dx} \cdot  \sqrt{{2\pi }/{20}}{\int }_{-\infty }^{\infty }f\left( {y; - {0.1},{20}^{-1}}\right) {dy}
\]

\[
= \sqrt{{2\pi }/{90}}\sqrt{{2\pi }/{120}} + {0.5}\sqrt{{2\pi }/{180}}\sqrt{{2\pi }/{20}}
\]

\[
\approx  {0.1128199}
\]

于是令

\[
g\left( {x,y}\right)  = \widetilde{g}\left( {x,y}\right) /{0.1128199}
\]

\[
= {0.5358984f}\left( {x; - {0.4},{90}^{-1}}\right) f\left( {y;{0.5},{120}^{-1}}\right)
\]

\[
+ {0.4641016f}\left( {x;{0.5},{180}^{-1}}\right) f\left( {y; - {0.1},{20}^{-1}}\right) \text{ , }
\]

\[
- \infty  < x < \infty , - \infty  < y < \infty ,
\]

用复合抽样法对 \(g\left( {x,y}\right)\) 抽样, \(N = {10000}\) 时一次模拟得到的 \({\widehat{I}}_{3} = {0.126}\) ,从这一次模拟估计的 \({\widehat{I}}_{3}\) 的渐近方差为 \({6.84} \times  {10}^{-8}\) ,重复模拟 \(B = {100}\) 次,则 100 个 \({\widehat{I}}_{3}\) 的平均值为 0.1258,样本方差为 \({5.37} \times  {10}^{-8}\) 。 \({\widehat{I}}_{2}\) 的样本方差是 \({\widehat{I}}_{3}\) 的样本方差的 363 倍,如果要达到相同的估计方差, 两种方法的样本量相差三百多倍。

例 3.2.5. 标准化的重要抽样法在贝叶斯统计推断中有重要作用。例如,设独立的观测样本 \({Y}_{j}\) 服从如下的贝塔一二项分布:

\[
f\left( {{y}_{j} \mid  K,\eta }\right)  = P\left( {{Y}_{j} = {y}_{j}}\right)  = \left( \begin{matrix} {n}_{j} \\  {y}_{j} \end{matrix}\right) \frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) },{y}_{j} = 0,1,\ldots ,{n}_{j}, \tag{3.26}
\]

其中 \(B\left( {\cdot , \cdot  }\right)\) 是贝塔函数, \({n}_{j}\) 为已知的正整数, \(K > 0,0 < \eta  < 1\) 为未知参数。贝塔一二项分布用于描述比二项分布更为分散的随机变量分布。按照贝叶斯统计的做法,假设参数 \(\left( {K,\eta }\right)\) 也是随机变量,具有所谓的 “先验分布”,假设 \(\left( {K,\eta }\right)\) 有如下的 “无信息” 先验分布密度:

\[
\pi \left( {K,\eta }\right)  \propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }, \tag{3.27}
\]

则 \(\left( {K,\eta }\right)\) 有如下的 “后验密度”:

\[
\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)  \propto  \pi \left( {K,\eta }\right) \mathop{\prod }\limits_{{j = 1}}^{n}f\left( {{y}_{j} \mid  K,\eta }\right)
\]

\[
\propto  \frac{1}{{\left( 1 + K\right) }^{2}}\frac{1}{\eta \left( {1 - \eta }\right) }\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {{K\eta } + {y}_{j},K\left( {1 - \eta }\right)  + {n}_{j} - {y}_{j}}\right) }{B\left( {{K\eta },K\left( {1 - \eta }\right) }\right) }. \tag{3.28}
\]

设要求 \(E\left( {\log K \mid  \mathbf{Y}}\right)  = {\int }_{0}^{\infty }\log K\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right) {dK}\) 的值。

如果可以从后验密度 \(\widetilde{p}\left( {K,\eta  \mid  \mathbf{Y}}\right)\) 直接抽样,可以用平均值法估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) ,但从 (3.28) 来看很难直接抽样。为此,使用标准化的重要抽样法。为了解除 \(\left( {K,\eta }\right)\) 的取值限制,作变换 \(\alpha  = \log K,\beta  = \log \frac{\eta }{1 - \eta }\) ,则 \(\alpha ,\beta  \in  \left( {-\infty ,\infty }\right)\) ,而(3.28)对应的 \(\left( {\alpha ,\beta }\right)\) 的后验密度为:

\[
p\left( {\alpha ,\beta  \mid  \mathbf{Y}}\right)  \propto  \frac{{e}^{\alpha }}{{\left( 1 + {e}^{\alpha }\right) }^{2}}\mathop{\prod }\limits_{{j = 1}}^{n}\frac{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }} + {y}_{j},\frac{{e}^{\alpha }}{1 + {e}^{\beta }} + {n}_{j} - {y}_{j}}\right) }{B\left( {\frac{{e}^{\alpha }}{1 + {e}^{-\beta }},\frac{{e}^{\alpha }}{1 + {e}^{\beta }}}\right) }. \tag{3.29}
\]

取值无限制的随机变量试抽样密度经常使用自由度较小的 \(\mathrm{t}\) 分布,比如 \(\mathrm{t}\left( 4\right)\) 分布,设 \(\mathrm{t}\left( 4\right)\) 分布密度函数为 \(g\left( \cdot \right)\) ,用独立的 \(\mathrm{t}\left( 4\right)\) 分布生成 \(\left( {\alpha ,\beta }\right)\) 的试抽样样本 \(\left( {{\alpha }_{i},{\beta }_{i}}\right) ,i = 1,2,\ldots ,N\) ,可以估计 \(E\left( {\log K \mid  \mathbf{Y}}\right)\) 为

\[
\widehat{\alpha } = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{\alpha }_{i}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}{\mathop{\sum }\limits_{{i = 1}}^{N}\frac{p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right) }{g\left( {\alpha }_{i}\right) g\left( {\beta }_{i}\right) }}. \tag{3.30}
\]

其中的 \(p\left( {{\alpha }_{i},{\beta }_{i} \mid  \mathbf{Y}}\right)\) 只要用(3.29)的右侧计算,因为分子和分母的归一化常数可以消掉。

\section*{3.2.5 分层抽样法}

用平均值法计算 \({\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) ,若 \(h\left( \mathbf{x}\right)\) 在 \(C\) 内取值变化范围大则估计方差较大。重要抽样法选取了与 \(f\left( x\right)\) 形状相似但是容易抽样的密度 \(g\left( x\right)\) 作为试投密度,大大提高了精度,但是这样的 \(g\left( \mathbf{x}\right)\) 有时难以找到。

如果把 \(C\) 上的积分分解为若干个子集上的积分,使得 \(h\left( \mathbf{x}\right)\) 在每个子集上变化不大,分别计算各个子集上的积分再求和,可以提高估计精度。这种方法与 \(§{2.2.5}\) 的复合抽样法类似, 叫做分层抽样法。这也是抽样调查中的重要技术。

例 3.2.6. 对函数

\[
h\left( x\right)  = \left\{  \begin{array}{ll} 1 + \frac{x}{10}, & 0 \leq  x \leq  {0.5} \\   - 1 + \frac{x}{10}, & {0.5} < x \leq  1 \end{array}\right.
\]

求定积分

\[
I = {\int }_{0}^{1}h\left( x\right) {dx}
\]

可以得 \(I\) 的精确值为 \(I = {0.05}\) 。我们用平均值法和分层抽样法来估计 \(I\) 并比较精度。

在 \(\left\lbrack  {0,1}\right\rbrack\) 区间随机抽取 \(N\) 点用平均值法得 \({\widehat{I}}_{2}\) ,其渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{\operatorname{Var}\left( {h\left( U\right) }\right) }{N} = \frac{143}{150N} \approx  \frac{0.9533}{N}.
\]

把 \(I\) 拆分为 \(\left\lbrack  {0,{0.5}}\right\rbrack\) 和 \(\left\lbrack  {{0.5},1}\right\rbrack\) 上的积分,即

\[
I = a + b = {\int }_{0}^{0.5}h\left( x\right) {dx} + {\int }_{0.5}^{1}h\left( x\right) {dx},
\]

对 \(a\) 和 \(b\) 分别用平均值法,得

\[
\widehat{a} = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = 1}}^{{N/2}}h\left( {{0.5}{U}_{i}}\right)  = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = 1}}^{{N/2}}\left( {1 + {0.05}{U}_{i}}\right) ,
\]

\[
\widehat{b} = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = \left( {N/2}\right)  + 1}}^{N}h\left( {{0.5} + {0.5}{U}_{i}}\right)  = \frac{0.5}{N/2}\mathop{\sum }\limits_{{i = \left( {N/2}\right)  + 1}}^{N}\left( {-1 + {0.05} + {0.05}{U}_{i}}\right) ,
\]

\[
{\widehat{I}}_{5} = \widehat{a} + \widehat{b},
\]

则分层抽样法结果 \({\widehat{I}}_{5}\) 的渐近方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{5}\right)  = \operatorname{Var}\left( {\widehat{a} + \widehat{b}}\right)  = \operatorname{Var}\left( \widehat{a}\right)  + \operatorname{Var}\left( \widehat{b}\right)
\]

\[
= {0.25}\frac{\operatorname{Var}\left( {1 + {0.05U}}\right) }{N/2} + {0.25}\frac{\operatorname{Var}\left( {-{0.95} + {0.05U}}\right) }{N/2} = \frac{1/{4800}}{N},
\]

分层后的估计方差远小于不分层的结果, 可以节省样本量约 4500 倍。

一般地,设积分 \(I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x}\) 可以分解为 \(m\) 个不交的子集 \({C}_{j}\) 上的积分,即

\[
I = {\int }_{C}h\left( \mathbf{x}\right) d\mathbf{x} = {\int }_{{C}_{1}}h\left( \mathbf{x}\right) d\mathbf{x} + {\int }_{{C}_{2}}h\left( \mathbf{x}\right) d\mathbf{x} + \cdots  + {\int }_{{C}_{m}}h\left( \mathbf{x}\right) d\mathbf{x}
\]

在 \({C}_{j}\) 投 \({n}_{j}\) 个随机点 \({X}_{ji} \sim  \mathrm{U}\left( {C}_{j}\right) ,i = 1,\ldots ,{n}_{j}\) ,则 \(I\) 的 \(m\) 个部分可以分别用平均值法估计,由此得 \(I\) 的分层估计为

\[
{\widehat{I}}_{5} = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{V\left( {C}_{j}\right) }{{n}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{n}_{j}}h\left( {X}_{ji}\right)
\]

记 \({\sigma }_{j}^{2} = \operatorname{Var}\left( {h\left( {X}_{j1}\right) }\right)\) ,划分子集时应使每一子集内 \(h\left( \cdot \right)\) 变化不大,即 \({\sigma }_{j}^{2}\) 较小。这时

\[
\operatorname{Var}\left( {\widehat{I}}_{5}\right)  = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{V}^{2}\left( {C}_{j}\right) {\sigma }_{j}^{2}}{{n}_{j}}
\]

若 \({\sigma }_{j}^{2}\) 可估计,应取 \({n}_{j}\) 使

\[
{n}_{j} \propto  V\left( {C}_{j}\right) {\sigma }_{j}, \tag{3.31}
\]

即

\[
{n}_{j} = N\frac{V\left( {C}_{j}\right) {\sigma }_{j}}{\mathop{\sum }\limits_{{k = 1}}^{m}V\left( {C}_{k}\right) {\sigma }_{k}},j = 1,2,\ldots ,m
\]

这样取的样本量 \(\left( {{n}_{1},{n}_{2},\ldots ,{n}_{m}}\right)\) 在所有满足 \({n}_{1} + {n}_{2} + \cdots  + {n}_{m} = N\) 的取法中使得渐近方差最小 (见习题1)。

在分层抽样法中, 划分了子集后, 每一子集上的积分也可用重要抽样法计算。

分层抽样法也可以用在求随机变量函数期望的问题中。设 \(X\) 为随机变量,要求 \(X\) 的函数 \(h\left( X\right)\) 的数学期望 \(\theta  = {Eh}\left( X\right)\) 。假设存在离散型随机变量 \(Y,{p}_{j} = P\left( {Y = {y}_{j}}\right) ,j = 1,2,\ldots ,m\) , 在 \(Y = {y}_{j}\) 条件下可以从 \(X\) 的条件分布抽样,则

\[
E\left\lbrack  {h\left( X\right) }\right\rbrack   = E\{ E\left\lbrack  {h\left( X\right)  \mid  Y}\right\rbrack  \}  = \mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  {p}_{j}, \tag{3.32}
\]

如果在 \(Y = {y}_{j}\) 条件下生成 \(X\) 的 \({N}_{j} = N{p}_{j}\) 个抽样值,设为 \({X}_{i}^{\left( j\right) },i = 1,2,\ldots ,{N}_{j}\) ,则可以用 \(\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right)\) 估计 \(E\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) ,估计 \(\theta\) 为

\[
\widehat{\theta } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right) {p}_{j} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}\mathop{\sum }\limits_{{i = 1}}^{{N{p}_{j}}}h\left( {X}_{i}^{\left( j\right) }\right) , \tag{3.33}
\]

这是 \(\theta\) 的无偏和强相合估计,且估计方差

\[
\operatorname{Var}\left( \widehat{\theta }\right)  = \frac{1}{{N}^{2}}\mathop{\sum }\limits_{{j = 1}}^{m}N{p}_{j}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack
\]

\[
= \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  {p}_{j} \tag{3.34}
\]

\[
= \frac{1}{N}E\{ \operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y}\right\rbrack  \}  \leq  \frac{1}{N}\operatorname{Var}\left\lbrack  {h\left( X\right) }\right\rbrack  , \tag{3.35}
\]

比直接用平均值法估计 \({Eh}\left( X\right)\) 的方差小。这里用到了条件方差的性质

\[
\operatorname{Var}\left( X\right)  = E\left\lbrack  {\operatorname{Var}\left( {X \mid  Y}\right) }\right\rbrack   + \operatorname{Var}\left\lbrack  {E\left( {X \mid  Y}\right) }\right\rbrack   \geq  E\left\lbrack  {\operatorname{Var}\left( {X \mid  Y}\right) }\right\rbrack  , \tag{3.36}
\]

如果 \(Y\) 与 \(X\) 独立则 \(E\left( {X \mid  Y}\right)  = {EX},\operatorname{Var}\left\lbrack  {E\left( {X \mid  Y}\right) }\right\rbrack   = 0\) ,这时分层抽样法比平均值法没有改进。从(3.34)可以看出,如果第 \(j\) 层样本的函数 \(h\left( {X}_{i}^{\left( j\right) }\right) ,i = 1,2,\ldots ,N{p}_{j}\) 的样本方差为 \({S}_{j}^{2}\) , 则 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 的一个无偏估计是

\[
\overset{⏜}{\operatorname{Var}\left( \widehat{\theta }\right) } = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{m}{S}_{j}^{2}{p}_{j} \tag{3.37}
\]

公式(3.33)取第 \(j\) 层样本数 \({N}_{j} = N{p}_{j}\) ,仅考虑了 \(Y\) 的取值分布,而未考虑 \(X \mid  Y = {y}_{j}\) 的条件分布情况。类似于(3.31)和(3.32),应该对 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) 较大的层取更多的样本。使得估计方差最小的分层样本量分配满足 \({N}_{j} \propto  {p}_{j}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }\) ,即

\[
{N}_{j} = N\frac{{p}_{j}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }}{\mathop{\sum }\limits_{{k = 1}}^{m}{p}_{k}\sqrt{\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{k}}\right\rbrack  }}. \tag{3.38}
\]

在 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) 未知的时候,可以预先抽取一个小的样本估计 \(\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack\) ,然后按估计的最优 \({N}_{j}\) 分配各层的样本量。采用(3.38)的分层样本量后,

\[
\widehat{\theta } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{1}{{N}_{j}}\mathop{\sum }\limits_{{i = 1}}^{{N}_{j}}h\left( {X}_{i}^{\left( j\right) }\right) {p}_{j}, \tag{3.39}
\]

\[
\operatorname{Var}\left( \widehat{\theta }\right)  = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{p}_{j}^{2}\operatorname{Var}\left\lbrack  {h\left( X\right)  \mid  Y = {y}_{j}}\right\rbrack  }{{N}_{j}}, \tag{3.40}
\]

于是 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 的估计为

\[
\widehat{\operatorname{Var}\left( \widehat{\theta }\right) } = \mathop{\sum }\limits_{{j = 1}}^{m}\frac{{p}_{j}^{2}{S}_{j}^{2}}{{N}_{j}}, \tag{3.41}
\]

其中 \({S}_{j}^{2}\) 是第 \(j\) 层样本函数 \(\left\{  {h\left( {X}_{i}^{\left( j\right) }\right) ,i = 1,2,\ldots ,{N}_{j}}\right\}\) 的样本方差。

分层抽样法的本质是把 \(X\) 的值相近的抽样分入一层,使得同层的 \(X\) 条件方差较小,从而减小估计方差。

例 3.2.7. 设 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,要估计 \(\theta  = {Eh}\left( U\right)  = {\int }_{0}^{1}h\left( x\right) {dx}\) 。令 \(Y = \operatorname{ceil}\left( {mU}\right)\) ,即当且仅当 \(\frac{j - 1}{m} < U \leq  \frac{j}{m}\) 时 \(Y = j,j = 1,2,\ldots ,m\) ,可以按照 \(Y\) 分层抽样估计 \(\theta\) :

\[
\theta  = E\left\lbrack  {h\left( U\right) }\right\rbrack   = \mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( U\right)  \mid  Y = j}\right\rbrack  P\left( {Y = j}\right)  \tag{3.42}
\]

\[
= \frac{1}{m}\mathop{\sum }\limits_{{j = 1}}^{m}E\left\lbrack  {h\left( U\right)  \mid  Y = j}\right\rbrack  , \tag{3.43}
\]

易见 \(Y = j\) 条件下 \(U\) 服从 \(\left( {\frac{j - 1}{m},\frac{j}{m}}\right)\) 上的均匀分布,设 \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 是 \(\mathrm{U}\left( {0,1}\right)\) 的独立抽样, 则用分层抽样法取每层 \({N}_{j} = 1\) 估计 \(\theta  = {Eh}\left( U\right)\) 为

\[
\widehat{\theta } = \frac{1}{m}\mathop{\sum }\limits_{{j = 1}}^{m}h\left( \frac{j - 1 + {U}_{j}}{m}\right) . \tag{3.44}
\]

\section*{3.3 方差缩减方法}

随机模拟方法虽然有着适用性广、方法简单的优点, 但是又有精度低、计算量大的缺点, 一整套模拟算几天几夜也是常有的事情。如果能成倍地减小随机模拟误差方差, 就可以有效地节省随机模拟时间, 有些情况下可以把耗时长到不具有可行性的模拟计算 (比如几个月) 缩短到可行 (比如几天)。

前一节的关于定积分计算的重要抽样法、分层抽样法都是降低随机模拟误差方差的重要方法, 也可以用在一般的模拟问题中。本节介绍一些其它的方差缩减技巧。我们以随机变量 \(X\) 的期望 \(\theta  = {EX}\) 的估计为例,目标是降低 \(\theta\) 的估计量的渐近方差。

\section*{3.3.1 控制变量法}

设要估计随机变量 \(X\) 的期望 \(\theta  = {EX}\) ,从 \(X\) 中抽取 \(N\) 个独立样本值 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) , 用样本平均值 \(\bar{X} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{X}_{i}\) 估计 \({EX}\) 。为了提高精度可以利用辅助信息。设有另外的随机变量 \(Y\) 满足

\[
{EY} = 0,\;\operatorname{Cov}\left( {X,Y}\right)  < 0
\]

令 \(Z = X + Y\) ,则

\[
E\left( Z\right)  = \theta ,\;\operatorname{Var}\left( Z\right)  = \operatorname{Var}\left( X\right)  + \operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right) ,
\]

只要 \(\operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right)  < 0\) 则 \(\operatorname{Var}\left( Z\right)  < \operatorname{Var}\left( X\right)\) ,如果有(X, Y)成对的抽样 \(\left( {{X}_{i},{Y}_{i}}\right) ,i =\)  \(1,2,\ldots ,n\) ,令 \({Z}_{i} = {X}_{i} + {Y}_{i}\) ,则用 \(\bar{Z}\) 来估计 \(\theta  = {EX} = {EZ}\) 的渐近方差就比用 \(\bar{X}\) 估计 \(I\) 的渐近方差减小了。

为了最好地利用 \(Y\) 与 \(X\) 的相关性,令

\[
Z\left( b\right)  = X + {bY},
\]

则

\[
{EZ}\left( b\right)  = {EX} = \theta ,
\]

\[
\operatorname{Var}\left( {Z\left( b\right) }\right)  = \operatorname{Var}\left( X\right)  + {2b}\operatorname{Cov}\left( {X,Y}\right)  + {b}^{2}\operatorname{Var}\left( Y\right) ,
\]

求 \(\operatorname{Var}\left( {Z\left( b\right) }\right)\) 关于 \(b\) 的最小值点,得

\[
b =  - \operatorname{Cov}\left( {X,Y}\right) /\operatorname{Var}\left( Y\right)  =  - {\rho }_{X,Y}\sqrt{\operatorname{Var}\left( X\right) /\operatorname{Var}\left( Y\right) },
\]

这时

\[
\operatorname{Var}\left( {Z\left( b\right) }\right)  = \left( {1 - {\rho }_{X,Y}^{2}}\right) \operatorname{Var}\left( X\right)  \leq  \operatorname{Var}\left( X\right) ,
\]

可见只要能找到零均值随机变量 \(Y\) 使得 \({\rho }_{X,Y} \neq  0\) 就可以减小 \({EX}\) 的估计方差。 \(Y\) 和 \(X\) 的相关性越强, 改善幅度越大。这种减小随机模拟误差方差的方法叫做控制变量法。

实际中 \({\rho }_{X,Y}\) 和 \(\operatorname{Var}\left( X\right) ,\operatorname{Var}\left( Y\right)\) 可能是未知的,可以先模拟一个小的样本估计 \({\rho }_{X,Y}\) 和 \(\operatorname{Var}\left( X\right) ,\operatorname{Var}\left( Y\right)\) 从而获得 \(b\) 的估计值。

控制变量法中要求控制变量 \(Y\) 与 \(X\) 相关且 \({EY} = 0\) 。如果 \({EY} \neq  0\) 但 \({EY} = {\mu }_{Y}\) 已知, 只要用 \(Y - {\mu }_{Y}\) 代替 \(Y\) ,这需要能预先知道 \(Y\) 的期望值的真值。另一种情况是 \({EY} = {EX}\) 未知, \(Y\) 与 \(X\) 相关,这时令 \(Z = {\alpha X} + \left( {1 - \alpha }\right) Y,\alpha  \in  \left\lbrack  {0,1}\right\rbrack\) ,则 \({EZ} = {EX} = \theta\) ,可以求 \(\alpha\) 使得 \(\operatorname{Var}\left( Z\right)\) 最小。容易知道当 \(\alpha  = \operatorname{Cov}\left( {Y,Y - X}\right) /\operatorname{Var}\left( {Y - X}\right)\) 时 \(\operatorname{Var}\left( Z\right)\) 最小。

例 3.3.1. 设要估计 \(I = {\int }_{0}^{1}{e}^{t}{dt}\) 。当然,可以得到积分真值为 \(e - 1\) ,这里用来演示控制变量法的优势。

设 \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {e}^{U}\) ,则 \(I = E{e}^{U} = {EX}\) ,可以用平均值法估计 \(\mathrm{I}\) 为

\[
{\widehat{I}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{e}^{{U}_{i}}. \tag{3.45}
\]

其方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{1}{N}\operatorname{Var}\left( {e}^{U}\right)  = \frac{1}{N}\left( {-\frac{1}{2}{e}^{2} + {2e} - \frac{3}{2}}\right)  \approx  \frac{0.2420}{N}. \tag{3.46}
\]

令 \(Y = U - \frac{1}{2}\) ,则 \({EY} = 0,X\) 与 \(Y\) 正相关,可以计算出 \(\operatorname{Cov}\left( {X,Y}\right)  \approx  {0.14086},\operatorname{Var}\left( Y\right)  = 1/{12}\) (更复杂的问题中可能需要从一个小的随机抽样中近似估计),于是 \(b =  - \operatorname{Cov}\left( {X,Y}\right) /\operatorname{Var}\left( Y\right)  =\) -1.690,对 \(Z\left( b\right)  = {e}^{U} - {1.690}\left( {U - \frac{1}{2}}\right)\) 有 \(\operatorname{Var}\left( {Z\left( b\right) }\right)  = \left\lbrack  {1 - {\rho }_{X,Y}^{2}}\right\rbrack  \operatorname{Var}\left( X\right)  = \left( {1 - {0.9919}^{2}}\right) \operatorname{Var}\left( X\right)  =\)  \({0.016}\operatorname{Var}\left( X\right)  = {0.0039}\) ,用控制变量法估计 \(I\) 为

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\left\lbrack  {{e}^{{U}_{i}} - {1.690}\left( {{U}_{i} - \frac{1}{2}}\right) }\right\rbrack  .
\]

\({\widehat{I}}_{1}\) 的方差比控制变量法 \({\widehat{I}}_{2}\) 的方差大 60 倍以上。

例 3.3.2. (系统可靠性估计) 考虑由 \(n\) 个部件组成的一个系统,用 \({S}_{i}\) 表示第 \(i\) 个部件是否正常工作,1 表示正常工作,0 表示失效。设 \({S}_{i} \sim  \mathrm{B}\left( {1,{p}_{i}}\right)\) 且各 \({S}_{i}\) 相互独立。用 \(Y\) 表示系统是否工作正常,1 表示工作正常,0 表示系统失效。设 \(Y\) 为 \({S}_{1},{S}_{2},\ldots ,{S}_{n}\) 的函数 \(\phi \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) 且 \(\phi\) 关于每个 \({S}_{i}\) 是单调不减,称 \(\phi\) 为系统的结构函数。令 \(R = P\left( {Y = 1}\right)  = {E\phi }\left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) , 称 \(R\) 为系统可靠度。

例如, \(\phi \left( {{s}_{1},{s}_{2},\ldots ,{s}_{n}}\right)  = \mathop{\prod }\limits_{{i = 1}}^{n}{s}_{i}\) ,则当且仅当所有部件正常工作时系统才正常工作,这样的系统称为串联系统, 这时系统可靠度为

\[
R = P\left( {{S}_{1} = 1,{S}_{2} = 1,\ldots ,{S}_{n} = 1}\right)  = \mathop{\prod }\limits_{{i = 1}}^{n}P\left( {{S}_{i} = 1}\right)  = {p}_{1}{p}_{2}\ldots {p}_{n}. \tag{3.47}
\]

在系统比较简单的情况下 (如串联、并联),可以给出用 \({p}_{1},{p}_{2},\ldots ,{p}_{n}\) 表示 \(R\) 的表达式。但是更复杂的系统则很难写出 \(R\) 的表达式,这时可以用随机模拟方法估计 \(R\) 。记 \(\mathbf{S} = \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right) ,\mathbf{X} = \phi \left( \mathbf{S}\right)\) ,对 \(\mathbf{S}\) 独立抽取 \(N\) 个点 \({\mathbf{S}}^{\left( j\right) } = \left( {{S}_{1}^{\left( j\right) },{S}_{2}^{\left( j\right) },\ldots ,{S}_{n}^{\left( j\right) }}\right) ,j =\)  \(1,2,\ldots ,N,R\) 可以用平均值法估计为

\[
{\widehat{R}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\phi \left( {\mathbf{S}}^{\left( j\right) }\right) . \tag{3.48}
\]

令 \(Y = \mathop{\sum }\limits_{{i = 1}}^{n}\left( {{S}_{i} - {p}_{i}}\right)\) ,则 \({EY} = 0,Y\) 与 \(X\) 正相关。用一个小的抽样先近似估计 \(\operatorname{Cov}\left( {X,Y}\right) ,\operatorname{Var}\left( Y\right)\) 得到 \(b\) 的近似值,可以得到方差缩减的估计量

\[
{\widehat{R}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\left\lbrack  {\phi \left( {\mathbf{S}}^{\left( j\right) }\right)  + b\mathop{\sum }\limits_{{i = 1}}^{n}\left( {{S}_{i}^{\left( j\right) } - {p}_{i}}\right) }\right\rbrack  .
\]

\section*{3.3.2 对立变量法}

控制变量法需要知道控制变量 \(Y\) 的期望值真值,并精确或近似知道 \(\operatorname{Cov}\left( {X,Y}\right)\) 和 \(\operatorname{Var}\left( Y\right)\) 。对立变量法的要求比较简单。

模拟中经常使用均匀分布随机数 \(U\) 变换产生的随机数 \(X = g\left( U\right)\) 。下面的定理说明,如果变换 \(g\left( \cdot \right)\) 是单调的,则随机变量 \(Y = g\left( {1 - U}\right)\) 就是与 \(X\) 负相关的。注意 \(g\left( {1 - U}\right)\) 与 \(g\left( U\right)\) 同分布所以 \({EY} = {EX}\) 。

定理 3.3.1. 设 \(g\) 为单调函数, \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,则 \(\operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) 。

证明. \(\forall {u}_{1},{u}_{2} \in  \left\lbrack  {0,1}\right\rbrack\) ,由 \(g\) 单调可知

\[
\left( {g\left( {u}_{1}\right)  - g\left( {u}_{2}\right) }\right) \left( {g\left( {1 - {u}_{1}}\right)  - g\left( {1 - {u}_{2}}\right) }\right)  \leq  0
\]

设 \({U}_{2}\) 服从 \(\mathrm{U}\left( {0,1}\right)\) 且与 \(U\) 独立,令 \({X}_{1} = g\left( U\right) ,{Y}_{1} = g\left( {1 - U}\right) ,{X}_{2} = g\left( {U}_{2}\right) ,{Y}_{2} = g\left( {1 - {U}_{2}}\right)\) , 则 \({X}_{1},{Y}_{1},{X}_{2},{Y}_{2}\) 的分布相同,且

\[
E\left( {{X}_{1} - {X}_{2}}\right) \left( {{Y}_{1} - {Y}_{2}}\right)
\]

\[
= \operatorname{Cov}\left( {{X}_{1} - {X}_{2},{Y}_{1} - {Y}_{2}}\right)
\]

\[
= \operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)  + \operatorname{Cov}\left( {{X}_{2},{Y}_{2}}\right)  - \operatorname{Cov}\left( {{X}_{1},{Y}_{2}}\right)  - \operatorname{Cov}\left( {{X}_{2},{Y}_{1}}\right)
\]

\[
= 2\operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)
\]

注意 \(\left( {{X}_{1} - {X}_{2}}\right) \left( {{Y}_{1} - {Y}_{2}}\right)  \leq  0\) 所以 \(\operatorname{Cov}\left( {{X}_{1},{Y}_{1}}\right)  \leq  0\) ,即 \(\operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) 。证毕。 定理3.3.1可以推广到如下情形。

定理 3.3.2. 设 \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 是关于每个自变量单调的函数, \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 相互独立,则 \(\operatorname{Cov}\left( {h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right) ,h\left( {1 - {U}_{1},1 - {U}_{2},\ldots ,1 - {U}_{n}}\right) }\right)  \leq  0\) 。

证明略去,参见 Ross \({\left( {2013}\right) }^{\left\lbrack  {34}\right\rbrack  }§{9.9}\) 。

对均匀随机数 \(U\) 最常见的变换是逆变换 \(X = {F}^{-1}\left( U\right)\) 。下面的定理给出了提高 \(I = {EX}\) 估计精度的方法。

定理 3.3.3 (对立变量法). 设 \(F\left( x\right)\) 为连续分布函数, \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {F}^{-1}\left( U\right) ,Y =\)  \({F}^{-1}\left( {1 - U}\right) ,Z = \frac{X + Y}{2}\) ,则 \(X\) 与 \(Y\) 同分布 \(F\left( x\right)\) 且 \(\operatorname{Cov}\left( {X,Y}\right)  \leq  0\) ,

\[
\operatorname{Var}\left( Z\right)  \leq  \frac{1}{2}\operatorname{Var}\left( X\right)
\]

证明. 因为 \(U\) 和 \(1 - U\) 同分布所以 \(X = {F}^{-1}\left( U\right)\) 和 \(Y = {F}^{-1}\left( {1 - U}\right)\) 同分布。由定理3.3.1, 令 \(g\left( \cdot \right)  = {F}^{-1}\left( \cdot \right)\) 可知 \(\operatorname{Cov}\left( {X,Y}\right)  = \operatorname{Cov}\left( {g\left( U\right) ,g\left( {1 - U}\right) }\right)  \leq  0\) ,从而

\[
\operatorname{Var}\left( Z\right)  = \frac{\operatorname{Var}\left( X\right)  + \operatorname{Var}\left( Y\right)  + 2\operatorname{Cov}\left( {X,Y}\right) }{4}
\]

\[
= \frac{\operatorname{Var}\left( X\right)  + \operatorname{Cov}\left( {X,Y}\right) }{2} \leq  \frac{1}{2}\operatorname{Var}\left( X\right)
\]

证毕。

根据定理 3.3.3 的结论,为了估计 \(I = {EX}\) ,产生 \({U}_{1},{U}_{2},\ldots ,{U}_{N}\) 后用

\[
{Z}_{i} = \frac{1}{2}\left( {{F}^{-1}\left( {U}_{i}\right)  + {F}^{-1}\left( {1 - {U}_{i}}\right) }\right) ,i = 1,2,\ldots ,N
\]

的样本平均值估计 \(I\) 可以提高精度,在不增加抽样个数的条件下把估计的随机误差方差降低到原来的 \(\frac{1}{2}\) 以下。这种提高随机模拟精度的方法叫做对立变量法。对立变量法不需要计算 \(X\) 和 \(Y\) 的方差及协方差的值,比控制变量法更简便易行。

例 3.3.3. 再次考虑例3.3.1的问题,估计 \(I = {\int }_{0}^{1}{e}^{t}{dt}\) 。下面用对立变量法改善原始的平均值法 \({\widehat{I}}_{1}\) 的估计方差。

设 \(U \sim  \mathrm{U}\left( {0,1}\right) ,X = {e}^{U}\) ,令 \(Y = {e}^{1 - U}\) ,用对立变量法估计 \(I\) 为

\[
{\widehat{I}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{{e}^{{U}_{i}} + {e}^{1 - {U}_{i}}}{2}, \tag{3.49}
\]

方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{3}\right)  = \frac{1}{N}\frac{\operatorname{Var}\left( {e}^{U}\right)  + \operatorname{Cov}\left( {{e}^{U},{e}^{1 - U}}\right) }{2} = \frac{1}{N}\left( {-\frac{3}{4}{e}^{2} + \frac{5}{2}e - \frac{5}{4}}\right)  \approx  \frac{0.003913}{N}, \tag{3.50}
\]

\({\widehat{I}}_{1}\) 的方差比对立变量法估计 \({\widehat{I}}_{3}\) 的方差大至少 60 倍,而 \({\widehat{I}}_{3}\) 的方差和例3.3.1中控制变量法估计量 \({\widehat{I}}_{2}\) 的方差相近。

对立变量法和控制变量法是类似做法, 一般不能结合使用。

例 3.3.4. 再次考虑例3.3.2的可靠度估计问题。用对立变量法改善估计方差。设 \(\left\{  {U}_{k}\right\}\) 为标准均匀分布随机数列, 取

\[
{S}_{i}^{\left( j\right) } = \left\{  \begin{array}{ll} 1 & \text{ 当 }{U}_{n\left( {j - 1}\right)  + i} \leq  {p}_{i}, \\  0 & \text{ 其它 } \end{array}\right.  \tag{3.51}
\]

则 \({S}_{i}^{\left( j\right) }\) 是 \({U}_{n\left( {j - 1}\right)  + i}\) 的单调非增函数。 \(R\) 用平均值法估计为

\[
{\widehat{R}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\phi \left( {{S}_{1}^{\left( j\right) },{S}_{2}^{\left( j\right) },\ldots ,{S}_{n}^{\left( j\right) }}\right) . \tag{3.52}
\]

利用对立变量法,令 \(h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right)  = \phi \left( {{S}_{1},{S}_{2},\ldots ,{S}_{n}}\right)\) ,则 \(h\) 关于每个自变量是单调非增函数, 于是

\[
\operatorname{Cov}\left( {h\left( {{U}_{1},{U}_{2},\ldots ,{U}_{n}}\right) ,h\left( {1 - {U}_{1},1 - {U}_{2},\ldots ,1 - {U}_{n}}\right) }\right)  \leq  0, \tag{3.53}
\]

估计系统可靠度 \(R\) 为

\[
{\widehat{R}}_{3} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}\frac{h\left( {{U}_{1}^{\left( j\right) },{U}_{2}^{\left( j\right) },\ldots ,{U}_{n}^{\left( j\right) }}\right)  + h\left( {1 - {U}_{1}^{\left( j\right) },1 - {U}_{2}^{\left( j\right) },\ldots ,1 - {U}_{n}^{\left( j\right) }}\right) }{2} \tag{3.54}
\]

就能比 \({\widehat{R}}_{1}\) 的误差方差至少降低 \(\frac{1}{2}\) 。

注意定理3.3.2中的 \({U}_{1},{U}_{2},\ldots ,{U}_{n}\) 仅要求独立,并未指定分布。事实上,此结果还可以推广。如果 \(\left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,n\) 独立,且 \({X}_{i}\) 和 \({Y}_{i}\) 负相关, \(h\left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 关于每个自变量单调不减,则 \(\operatorname{Cov}\left( {h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right) ,h\left( {{Y}_{1},{Y}_{2},\ldots ,{Y}_{n}}\right) }\right)  \leq  0\) ,可以类似地用控制变量法或对立变量法的做法构造方差缩减估计量。例如,设 \({X}_{i} \sim  \mathrm{N}\left( {{\mu }_{i},{\sigma }_{i}^{2}}\right)\) 相互独立,要估计 \(\theta  =\)  \({Eh}\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) ,其中 \(h\) 关于每个自变量单调不减,则 \(h\left( {2{\mu }_{1} - {X}_{1},2{\mu }_{2} - {X}_{2},\ldots ,2{\mu }_{n} - {X}_{n}}\right)\) 与 \(h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) 同分布,对

\[
Z = \frac{h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)  + h\left( {2{\mu }_{1} - {X}_{1},2{\mu }_{2} - {X}_{2},\ldots ,2{\mu }_{n} - {X}_{n}}\right) }{2}
\]

抽样用平均值估计 \(\theta\) 可以比仅对 \(h\left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) 抽样得到的估计量的方差缩减一半以上。

\section*{3.3.3 条件期望法}

进行统计估计时, 如果有额外的相关信息, 利用这样的信息可以提高估计精度。比如, 对随机变量 \(Y\) ,如果 \(Y\) 服从某种模型,在估计 \(I = {EY}\) 时应当尽量利用模型信息。

设变量 \(X\) 与 \(Y\) 不独立,根据 Rao-Blackwell 不等式:

\[
\operatorname{Var}\{ \mathrm{E}\left( {Y \mid  X}\right) \}  \leq  \operatorname{Var}\left( Y\right)
\]

又

\[
\mathrm{E}\{ \mathrm{E}\left( {Y \mid  X}\right) \}  = \mathrm{E}Y = I
\]

所以,对 \(Z = \mathrm{E}\left( {Y \mid  X}\right)\) 抽样,用 \(Z\) 的样本平均值来估计 \(I = \mathrm{E}Y\) 比直接用 \(Y\) 的样本平均值的精度更高。这种改善随机模拟估计精度的方法叫做条件期望法, 或 Rao-Blackwell 方法。

例 3.3.5. 设 \(X \sim  {p}_{1}\left( x\right) ,\varepsilon  \sim  \mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 且与 \(X\) 独立,

\[
Y = \psi \left( X\right)  + \varepsilon ,
\]

要估计 \(I = \mathrm{E}Y\) ,可以用条件分布抽样法对二元随机向量 \(\mathbf{Z} = \left( {X,Y}\right)\) 抽样产生 \(Y\) 的样本。从 \(p\left( \cdot \right)\) 抽样得 \({X}_{1},{X}_{2},\ldots ,{X}_{N}\) ,独立地从 \(\mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 抽样得 \({\varepsilon }_{1},{\varepsilon }_{2},\ldots ,{\varepsilon }_{N}\) ,令

\[
{Y}_{i} = \psi \left( {X}_{i}\right)  + {\varepsilon }_{i},i = 1,2,\ldots ,N
\]

然后用 \({Y}_{1},{Y}_{2},\ldots ,{Y}_{N}\) 的样本平均值估计 \({EY}\) :

\[
{\widehat{I}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{Y}_{i}
\]

估计方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{1}\right)  = \frac{\operatorname{Var}\left( {Y}_{1}\right) }{N} = \frac{\operatorname{Var}\left( {\psi \left( {X}_{1}\right) }\right) }{N} + \frac{{\sigma }^{2}}{N}.
\]

另一方面,注意 \(E\left( {Y \mid  X}\right)  = \psi \left( X\right)\) ,也可以只对 \(X\) 抽样然后用条件期望法估计 \({EY}\) :

\[
{\widehat{I}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{i}\psi \left( {X}_{i}\right) ,
\]

估计方差为

\[
\operatorname{Var}\left( {\widehat{I}}_{2}\right)  = \frac{\operatorname{Var}\left( {\psi \left( {X}_{1}\right) }\right) }{N} < \operatorname{Var}\left( {\widehat{I}}_{1}\right) .
\]

这个例子演示了条件期望法可以缩减误差方差的原因: 对 \(Y\) 抽样分成两步进行: 第一步对 \(X\) 抽样,第二步利用 \(Y \mid  X\) 的条件分布对 \(Y\) 抽样。所以使用 \(Y\) 的样本估计 \({EY}\) 包含了第二步对 \(Y\) 抽样的随机误差,而使用 \(X\) 的函数 \(E\left( {Y \mid  X}\right)  = \psi \left( X\right)\) 的抽样来估计 \({EY}\) 则避免了第二步对 \(Y\) 抽样引起的随机误差。

例 3.3.6. 继续考虑例3.1.1中用随机模拟方法估计 \(\pi\) 的问题。设(X, Y)服从正方形 \(D =\)  \({\left\lbrack  -1,1\right\rbrack  }^{2}\) 上的均匀分布,令 \(\eta  = 1\) 表示(X, Y)落入单位圆 \(C = \left\{  {\left( {x,y}\right)  : {x}^{2} + {y}^{2} \leq  1}\right\}  ,\eta  = 0\) 表示未落入单位圆,则 \(\eta  \sim  \mathrm{B}\left( {1,\pi /4}\right)\) 。设 \(\left( {{X}_{i},{Y}_{i}}\right) ,i = 1,2,\ldots ,N\) 是(X, Y)的独立重复抽样, \({\eta }_{i}\) 表示 \({X}_{i}^{2} + {Y}_{i}^{2} \leq  1\) 是否发生,则 \(\pi\) 的估计为

\[
{\widehat{\pi }}_{1} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}{\eta }_{i}
\]

方差为 \(\operatorname{Var}\left( {\widehat{\pi }}_{1}\right)  = \pi \left( {4 - \pi }\right) /N \approx  \frac{2.6968}{N}\) 。

用 \(\zeta  = E\left( {\eta  \mid  X}\right)\) 的样本代替 \(\eta\) 来估计 \({E\eta } = \pi /4\) ,则由

\[
E\left( {\eta  \mid  X = x}\right)  = P\left( {{X}^{2} + {Y}^{2} \leq  1 \mid  X = x}\right)  = P\left( {{Y}^{2} \leq  1 - {x}^{2}}\right)
\]

\[
= P\left( {-\sqrt{1 - {x}^{2}} \leq  Y \leq  \sqrt{1 - {x}^{2}}}\right)
\]

\[
= \sqrt{1 - {x}^{2}}\;\text{ (注意 }Y \sim  \mathrm{U}\left( {-1,1}\right) \text{ ) }
\]

可知 \(\zeta  = \sqrt{1 - {X}^{2}}\) 。其方差为

\[
\operatorname{Var}\left( \zeta \right)  = E\left( {1 - {X}^{2}}\right)  - {\left( E\zeta \right) }^{2} = \frac{2}{3} - \frac{{\pi }^{2}}{16},
\]

估计 \(\pi\) 为

\[
{\widehat{\pi }}_{2} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\sqrt{1 - {X}_{i}^{2}},
\]

方差为 \(\operatorname{Var}\left( {\widehat{\pi }}_{2}\right)  \approx  {0.7971}/N,\operatorname{Var}\left( {\widehat{\pi }}_{1}\right) /\operatorname{Var}\left( {\widehat{\pi }}_{2}\right)  \approx  {3.4}\) 。

另外, \(\zeta\) 仅依赖于 \({X}^{2}\) ,容易发现若 \(U \sim  \mathrm{U}\left( {0,1}\right)\) 则 \({X}^{2}\) 和 \({U}^{2}\) 同分布,所以可取 \(\xi  = \sqrt{1 - {U}^{2}}\) ,其中 \(U \sim  \mathrm{U}\left( {0,1}\right)\) 。这时,函数 \(h\left( u\right)  = \sqrt{1 - {u}^{2}}\) 是 \(u \in  \left( {0,1}\right)\) 的单调函数,可以利用对立变量法,构造 \(\pi\) 的估计量为

\[
{\widehat{\pi }}_{3} = \frac{4}{N}\mathop{\sum }\limits_{{i = 1}}^{N}\frac{\sqrt{1 - {U}_{i}^{2}} + \sqrt{1 - {\left( 1 - {U}_{i}\right) }^{2}}}{2},
\]

其中 \({U}_{1},{U}_{2},\ldots ,{U}_{N}\) 为 \(U\left( {0,1}\right)\) 随机数,则 \(\operatorname{Var}\left( {\widehat{\pi }}_{3}\right)  \approx  {0.11}/N,\operatorname{Var}\left( {\widehat{\pi }}_{1}\right) /\operatorname{Var}\left( {\widehat{\pi }}_{3}\right)  \approx  {25}\) 。

也可以用对立变量法来改进 \({\widehat{\pi }}_{2}\) 。令 \(W = {U}^{2} - \frac{1}{3}\) ,则 \({EW} = 0\) 且 \(W\) 与 \(\zeta\) 负相关。可以先进行一个小规模的模拟估计 \(\operatorname{Cov}\left( {\zeta ,W}\right)\) 和 \(\operatorname{Var}\left( W\right)\) 得到 \(b =  - \operatorname{Cov}\left( {\zeta ,W}\right) /\operatorname{Var}\left( W\right)\) 的近似值,用 \(\zeta \left( b\right)  = \zeta  + {bW}\) 代替 \(\zeta\) 进行抽样,可以减小 \({\widehat{\pi }}_{2}\) 的方差。

例 3.3.7. 设随机变量 \(\mathbf{Y} \sim  f\left( \mathbf{y}\right)\) ,为了估计 \(\theta  = {Eh}\left( \mathbf{Y}\right)\) ,经常使用标准化重要抽样法

\[
{\widehat{\theta }}_{1} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}h\left( {\mathbf{X}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{W}_{i}}, \tag{3.55}
\]

其中 \({\mathbf{X}}_{i}\) 为试抽样密度 \(g\left( \mathbf{x}\right)\) 的样本,权重 \({W}_{i} = f\left( {\mathbf{X}}_{i}\right) /g\left( {\mathbf{X}}_{i}\right)\) 。如果 \(\mathbf{x} = \left( {\mathbf{u},\mathbf{v}}\right) ,h\left( \mathbf{x}\right)  =\)  \({h}_{1}\left( \mathbf{u}\right)\) ,则只要对 \(\mathbf{X} = \left( {\mathbf{U},\mathbf{V}}\right)\) 的分量 \(\mathbf{U}\) 抽样得到 \({\mathbf{U}}_{i},i = 1,\ldots ,N\) 及新的权重 \({\widetilde{W}}_{i} =\)  \({f}_{\mathbf{U}}\left( {\mathbf{U}}_{i}\right) /{g}_{\mathbf{U}}\left( {\mathbf{U}}_{i}\right)\) ,其中 \({f}_{\mathbf{U}}\left( \mathbf{u}\right)  = \int f\left( {\mathbf{u},\mathbf{v}}\right) d\mathbf{v},{g}_{\mathbf{U}}\left( \mathbf{u}\right)  = \int g\left( {\mathbf{u},\mathbf{v}}\right) d\mathbf{v}\) ,则可估计 \(\theta\) 为

\[
{\widehat{\theta }}_{2} = \frac{\mathop{\sum }\limits_{{i = 1}}^{N}{\widetilde{W}}_{i}{h}_{1}\left( {\mathbf{U}}_{i}\right) }{\mathop{\sum }\limits_{{i = 1}}^{N}{\widetilde{W}}_{i}}, \tag{3.56}
\]

这时

\[
\operatorname{Var}\left( {\widetilde{W}}_{i}\right)  \leq  \operatorname{Var}\left( {W}_{i}\right) ,i = 1,2,\ldots ,N. \tag{3.57}
\]

事实上,

\[
\frac{{f}_{\mathbf{U}}\left( \mathbf{u}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) } = \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) }{dv} = \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{{g}_{\mathbf{U}}\left( \mathbf{u}\right) {g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) }{g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) d\mathbf{v}
\]

\[
= \int \frac{f\left( {\mathbf{u},\mathbf{v}}\right) }{g\left( {\mathbf{u},\mathbf{v}}\right) }{g}_{\mathbf{V} \mid  \mathbf{U}}\left( {\mathbf{v} \mid  \mathbf{u}}\right) d\mathbf{v} = {E}_{g}\left\{  {\left. \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right| \;\mathbf{U} = \mathbf{u}}\right\}   \tag{3.58}
\]

于是

\[
{\operatorname{Var}}_{g}\left\{  \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right\}   \geq  {\operatorname{Var}}_{g}\left\{  {{E}_{g}\left\lbrack  {\left. \frac{f\left( {\mathbf{U},\mathbf{V}}\right) }{g\left( {\mathbf{U},\mathbf{V}}\right) }\right| \;\mathbf{U}}\right\rbrack  }\right\}   = {\operatorname{Var}}_{g}\left\{  \frac{{f}_{\mathbf{U}}\left( \mathbf{U}\right) }{{g}_{\mathbf{U}}\left( \mathbf{U}\right) }\right\}  . \tag{3.59}
\]

\section*{3.3.4 随机数复用}

在统计研究中, 经常需要比较若干种统计方法的性能, 如偏差、方差、覆盖率等。除了努力获取理论结果以外,可以用随机模拟方法进行比较: 重复生成 \(N\) 组随机样本,对每个样本同时用不同统计方法计算结果,最后从 \(N\) 组结果比较不同方法的性能。这样比较时,并没有对每种方法单独生成 \(N\) 组样本,而是每个样本同时应用所有要比较的方法。这样不仅减少了计算量,而且在比较时具有更高的精度。

例 3.3.8. 对正态分布总体 \(X \sim  N\left( {\mu ,{\sigma }^{2}}\right)\) ,如果有样本 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) ,估计 \({\sigma }^{2}\) 有两种不同的公式:

\[
{\widehat{\sigma }}_{1}^{2} = \frac{1}{n - 1}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i} - \bar{X}\right) }^{2},{\widehat{\sigma }}_{2}^{2} = \frac{1}{n}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i} - \bar{X}\right) }^{2}. \tag{3.60}
\]

希望比较两个估计量的偏差和均方误差:

\[
{b}_{1} = E{\widehat{\sigma }}_{1}^{2} - {\sigma }^{2},\;{b}_{2} = E{\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}, \tag{3.61}
\]

\[
{s}_{1} = E{\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2},\;{s}_{2} = E{\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}. \tag{3.62}
\]

当然, 这个问题很简单, 可以得到偏差和均方误差的理论值:

\[
{b}_{1} = 0,\;{b}_{2} =  - \frac{1}{n}{\sigma }^{2}, \tag{3.63}
\]

\[
{s}_{1} = \frac{2{\sigma }^{4}}{n - 1},\;{s}_{2} = \frac{\left( {{2n} - 1}\right) {\sigma }^{4}}{{n}^{2}},\;{s}_{1} - {s}_{2} = \frac{\left( {{3n} - 1}\right) {\sigma }^{4}}{{n}^{2}\left( {n - 1}\right) }. \tag{3.64}
\]

我们用随机模拟来作比较。重复地生成 \(N\) 组样本 \(\left( {{X}_{1}^{\left( j\right) },{X}_{2}^{\left( j\right) },\ldots ,{X}_{n}^{\left( j\right) }}\right) ,j = 1,2,\ldots ,N\) 。对每组样本分别计算 \({\widehat{\sigma }}_{1j}^{2} = \frac{1}{n - 1}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i}^{\left( j\right) } - {\bar{X}}^{\left( j\right) }\right) }^{2}\) 和 \({\widehat{\sigma }}_{2j}^{2} = \frac{1}{n}\mathop{\sum }\limits_{{i = 1}}^{n}{\left( {X}_{i}^{\left( j\right) } - {\bar{X}}^{\left( j\right) }\right) }^{2}\) ,得到偏差和均方误差的估计

\[
{\widehat{b}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\widehat{\sigma }}_{1j}^{2} - {\sigma }^{2},\;{\widehat{b}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\widehat{\sigma }}_{2j}^{2} - {\sigma }^{2}, \tag{3.65}
\]

\[
{\widehat{s}}_{1} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\left( {\widehat{\sigma }}_{1j}^{2} - {\sigma }^{2}\right) }^{2},\;{\widehat{s}}_{2} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{\left( {\widehat{\sigma }}_{2j}^{2} - {\sigma }^{2}\right) }^{2}. \tag{3.66}
\]

容易看出, 两种方法使用相同的模拟样本得到的偏差、均方误差的估计精度与每种方法单独生成模拟样本得到的估计精度相同。

但是,如果要估计 \({\Delta s} = {s}_{1} - {s}_{2}\) ,利用相同的样本的估计精度更好。一般地,

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \operatorname{Var}\left( {\widehat{s}}_{1}\right)  + \operatorname{Var}\left( {\widehat{s}}_{2}\right)  - 2\operatorname{Cov}\left( {{\widehat{s}}_{1},{\widehat{s}}_{2}}\right) . \tag{3.67}
\]

如果每种方法使用不同的样本, 则 (3.67) 变成

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \frac{1}{N}\left\lbrack  {\operatorname{Var}\left( {\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\right)  + \operatorname{Var}\left( {\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\right) }\right\rbrack  . \tag{3.68}
\]

如果两种方法利用相同的样本计算, 则 (3.67) 变成

\[
\operatorname{Var}\left( {{\widehat{s}}_{1} - {\widehat{s}}_{2}}\right)  = \frac{1}{N}\left\lbrack  {\operatorname{Var}\left( {\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\right)  + \operatorname{Var}\left( {\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\right)  - 2\operatorname{Cov}\left( {{\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2},{\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}}\right) }\right\rbrack  .
\]

(3.69)

而 \({\left( {\widehat{\sigma }}_{1}^{2} - {\sigma }^{2}\right) }^{2}\) 与 \({\left( {\widehat{\sigma }}_{2}^{2} - {\sigma }^{2}\right) }^{2}\) 明显具有强正相关 (见习题10),所以两种方法针对相同样本计算时 \({\widehat{s}}_{1} - {\widehat{s}}_{2}\) 的方差比两种方法使用单独样本时的方差要小得多。这说明重复利用相同的随机数或样本往往可以提高比较的精度。

\section*{3.4 随机服务系统模拟 *}

我们在 \(§{3.1}\) 中讲到,较为复杂的随机模型往往难以进行彻底的理论分析,这时常常使用随机模拟方法产生模型的大量数据, 从产生的数据对模型进行统计推断。随机服务系统就是这样的一种模型, 经常需要利用随机模拟方法进行研究。

随机服务系统在我们日常生活、工业生产、科学技术、军事领域中是经常遇到的随机模型,比如,研究银行、理发店、商店、海关通道、高速路收费口等服务人员个数的设置和排队规则,研究计算机网络网关、移动网络的调度规则,等等。

在概率统计理论中排队论用来研究随机服务系统的数学模型, 可以用来设计适当的服务机构数目和排队规则。如下面的 \(\mathrm{M}/\mathrm{M}/1\) 排队系统。

例 3.4.1. 设某银行仅有一个柜员,并简单假设银行不休息。顾客到来间隔的时间服从独立的指数分布 \(\operatorname{Exp}\left( \lambda \right) (1/\lambda\) 为间隔时间的期望值),如果柜员正在为先前的顾客服务,新到顾客就排队等待,柜员为顾客服务的时间服从均值为 \(1/\mu\) 的指数分布,设 \(u = \lambda /\mu  < 1\) 。设 \({X}_{t}\) 表示 \(t\) 时刻在银行内的顾客数 (包括正在服务的和正在排队的),则 \({X}_{t}\) 是一个连续时马氏链。

这是一个生灭过程马氏链, 有理论结果当系统处于稳定状态时

\[
P\left( {{X}_{t} = i}\right)  = {u}^{i}\left( {1 - u}\right) ,i = 0,1,2,\ldots
\]

设随机变量 \(N\) 服从 \({X}_{t}\) 的平稳分布。于是银行中平均顾客数为

\[
{EN} = \frac{u}{1 - u}
\]

平均队列长度 \({EQ}\) 等于 \({EN}\) 减去平均正在服务人数,正在服务人数 \({Y}_{t}\) 为

\[
{Y}_{t} = \left\{  \begin{array}{ll} 1, & \text{ 当 }{X}_{t} > 0, \\  0, & \text{ 当 }{X}_{t} = 0 \end{array}\right.
\]

所以 \(E{Y}_{t} = P\left( {{Y}_{t} = 1}\right)  = 1 - P\left( {N = 0}\right)  = u\) ,于是平均队列长度为

\[
{EQ} = {EN} - E{Y}_{t} = \frac{u}{1 - u} - u = \frac{{u}^{2}}{1 - u},
\]

设顾客平均滞留时间为 \({ER}\) ,由关系式

\[
{EN} = \lambda  \cdot  {ER}
\]

可知平均滞留时间为

\[
{ER} = \frac{u}{\lambda \left( {1 - u}\right) } = \frac{1}{\mu  - \lambda },
\]

进一步分析还可以知道顾客滞留时间 \(R\) 服从均值为 \(1/\left( {\mu  - \lambda }\right)\) 的指数分布。

从上面的例子可以发现, 一个随机服务系统的模型应该包括如下三个要素:

\begin{itemize}
\item 输入过程: 比如, 银行的顾客到来的规律。
\end{itemize}

\begin{itemize}
\item 排队规则: 比如, 银行有多个柜员时, 顾客是选最短一队, 还是随机选一队, 还是统一排队等候叫号, 顾客等得时间过长后是否会以一定概率放弃排队。
\end{itemize}

\begin{itemize}
\item 服务机构:有多少个柜员,服务时间的分布等。
\end{itemize}

虽然某些随机服务系统可以进行严格理论分析得到各种问题的理论解, 但是, 随机服务系统中存在大量随机因素, 使得理论分析变得很困难以至于不可能。例如, 即使是上面的银行服务问题, 可能的变化因素就包括: 顾客到来用齐次泊松过程还是非齐次泊松过程, 柜员有多少个, 是否不同时间段柜员个数有变化, 柜员服务时间服从什么样的分布, 顾客排队按照什么规则,是否 VIP 顾客提前服务,顾客等候过长时会不会放弃排队,等等。包含了这么多复杂因素的随机服务系统的理论分析会变得异常复杂, 完全靠理论分析无法解决问题, 这时, 可以用随机模拟方法给出答案。

在模拟随机服务系统时, 可以按时间顺序记录发生的事件, 如顾客到来、顾客接受服务、 顾客结束服务等,这样的系统的模拟也叫做离散事件模拟。

离散事件模拟算法可以分为三类:

\begin{itemize}
\item 活动模拟, 把连续时间离散化为很小的单位, 比如平均几秒发生一个新事件时可以把时间精确到毫秒, 然后时钟每隔一个小时间单位前进一步并查看所有活动并据此更新系统状态。缺点是速度太慢。
\end{itemize}

\begin{itemize}
\item 事件模拟。仅在事件发生时更新时钟和待发生的事件集合。这样的方法不受编程语言功能的限制, 运行速度很快, 也比较灵活, 可以实现复杂逻辑, 但是需要自己管理的数据结构和逻辑结构比较复杂,算法编制相对较难。
\end{itemize}

\begin{itemize}
\item 过程模拟。把将要到来的各种事件看成是不同的过程, 让不同的过程并行地发生, 过程之间可以交换消息, 并在特殊的软件包或编程语言支持下自动更新系统状态。在计算机实现中需要借助于线程或与线程相似的程序功能。这类软件包有 \(\mathrm{C} +  +\) 语言软件包 \(\mathrm{C} +  +\) SIM 和 Python 语言软件包 SimPy。优点是系统逻辑的编码很直观,程序模块化, 需要用户自己管理的数据结构和逻辑结构少。过程模拟是现在更受欢迎的离散事件模拟方式。
\end{itemize}

事件模拟算法必须考虑的变量包括:

\begin{itemize}
\item 当前时刻 \(t\) ;
\end{itemize}

\begin{itemize}
\item 随时间而变化的计数变量,如 \(t\) 时刻时到来顾客人数、已离开人数;
\end{itemize}

\begin{itemize}
\item 系统状态, 比如是否有顾客正在接受服务、排队人数、队列中顾客序号。
\end{itemize}

这些变量仅在有事件发生时 (如顾客到来、顾客离开) 才需要记录并更新。为了持续模拟后续事件, 需要维护一个将要发生的事件的列表 (下一个到来时刻、下一个离开时刻), 列表仅需要在事件发生时进行更新。其它变量可以从这三种变量中推算出来,比如,顾客 \(i\) 在时间 \({t}_{1}\) 时到达并排队,在时间 \({t}_{2}\) 时开始接受服务,在时间 \({t}_{4}\) 时结束服务离开,则顾客 \(i\) 的滞留时间为 \({t}_{4} - {t}_{1}\) 。

例 3.4.2. 用事件模拟的方法来模拟例3.4.1。目的是估计平均滞留时间 \({ER}\) 。想法是,模拟生成各种事件发生的时间,模拟很长时间,丢弃开始的一段时间 \({T}_{0}\) 后,用 \({T}_{0}\) 后到达的顾客的总滞留时间除以 \({T}_{0}\) 后到达的顾客人数来估计平均滞留时间。设 \({T}_{0}\) 时间后每位顾客滞留时间的模拟值为 \({R}_{i},i = 1,2,\ldots ,m\) ,可以用 \(\left\{  {R}_{i}\right\}\) 作为随机变量 \(R\) 的样本来检验 \(R\) 的分布是否指数分布。

用事件模拟方法进行离散事件模拟的算法关键在于计算系统状态改变的时间,即各个事件的发生时间, 这个例子中就是顾客到来、顾客开始接受服务、顾客离开这样三种事件, 由此还可以得到每个顾客排队的时间和服务的时间。

在没有明确算法构思时, 可以从时间 0 开始在纸上按照模型规定人为地生成一些事件并人为地找到下一事件。这样可以找到要更新的数据结构和更新的程序逻辑。

下面的算法保持了一个将要发生的事件的集合, 在每次事件发生时更新时钟, 更新时钟时从事件集合中找到最早发生的事件进行处理, 并生成下一事件到事件集合中, 如此重复直到需要模拟的时间长度。

\HRule

\{初始化当前时钟 \(t \leftarrow  0\) ,柜员忙标志 \(B \leftarrow  0\) ,当前排队人数 \(L \leftarrow  0\) ,

最新到来顾客序号 \(i \leftarrow  0\) ,正在服务顾客序号 \(j \leftarrow  0\) ,已服务顾客数 \(n \leftarrow  0\}\)

从 \(\operatorname{Exp}\left( \lambda \right)\) 抽取 \(X\) ,设置下一顾客来到时间 \(A \leftarrow  X\)

repeat \{

\hspace*{1em} if \(\left( {B = 0\text{ or }\left( {B = 1\text{ and }A < E}\right) }\right) \left\{  {\# E}\right.\) 是正在服务的顾客结束时刻

\hspace*{2em} \(t \leftarrow  A\)

\hspace*{1em} \} else \{

\hspace*{2em} \(t \leftarrow  E\)

\hspace*{1em} \}

\hspace*{1em} if \(\left( {t > {T}_{1}}\right)\) break \(\# {T}_{1}\) 是预先确定的模拟时长

\hspace*{1em} if \(\left( {t =  = A}\right) \{ \#\) 待处理到达事件

\hspace*{2em} \(L \leftarrow  L + 1\)

\hspace*{2em} \(i \leftarrow  i + 1\) ,记录第 \(i\) 位顾客到来时间 \({a}_{i} \leftarrow  t\)

\hspace*{2em} 从 \(\operatorname{Exp}\left( \lambda \right)\) 抽取 \(X,A \leftarrow  t + X\)

\hspace*{2em} if \(\left( {B = 0}\right) \{ \#\) 不用排队,直接服务

\hspace*{3em} \(B \leftarrow  1,L \leftarrow  L - 1\)

\hspace*{3em} \(j \leftarrow  j + 1\) ,置第 \(j\) 位顾客开始服务时间 \({s}_{j} \leftarrow  t\)

\hspace*{3em} 从 \(\operatorname{Exp}\left( \mu \right)\) 抽取 \(Y\) ,置 \(E \leftarrow  t + Y\)

\hspace*{2em} \}

\hspace*{1em} \} else \{ \# 待处理结束服务事件

\hspace*{2em} \(B \leftarrow  0\)

\hspace*{2em} \(n \leftarrow  n + 1\) ,记录第 \(n\) 个顾客结束服务时间 \({e}_{n} \leftarrow  t\)

\hspace*{2em} if \(\left( {L > 0}\right) \{ \#\) 排队顾客开始服务

\hspace*{3em} \(L \leftarrow  L - 1\)

\hspace*{3em} \(B \leftarrow  1\)

\hspace*{3em} \(j \leftarrow  j + 1,{s}_{j} \leftarrow  t\)

\hspace*{2em} 从 \(\operatorname{Exp}\left( \mu \right)\) 抽取 \(Y\) ,置 \(E \leftarrow  t + Y\)

\hspace*{1em} \}

\}

\HRule

\}

\{令 \(I = \left\{  {i : {T}_{0} \leq  {s}_{i} \leq  {T}_{1}}\right\}\) ,求 \(\left\{  {{e}_{i} - {a}_{i},i \in  I}\right\}\) 的平均值作为 \({ER}\) 估计 \(\}\)

从这个算法可以看出, 事件模拟方法需要用户自己管理待处理事件集合与时钟, 算法设计难度较大。

用随机服务系统进行建模和模拟研究的一般步骤如下:

\begin{itemize}
\item 提出问题。比如, 银行中顾客平均滞留时间与相关参数的关系。
\end{itemize}

\begin{itemize}
\item 建立模型。比如, 设顾客到来服从泊松过程, 设每个顾客服务时间服从独立的指数分布。
\end{itemize}

\begin{itemize}
\item 数据收集与处理。比如,估计银行顾客到来速率,每个顾客平均服务时间,等等。
\end{itemize}

\begin{itemize}
\item 建立模拟程序。一般使用专用的模拟软件包或专用模拟语言编程。
\end{itemize}

\begin{itemize}
\item 模拟模型的正确性确认。
\end{itemize}

\begin{itemize}
\item 模拟试验和模拟结果分析。
\end{itemize}

离散事件模拟问题一般都比较复杂, 即使借助于专用模拟软件或软件包, 也很难确保实现的模拟算法与问题的实际设定是一致的。为此, 需要遵循一些提高算法可信度的一般规则。 算法程序一定要仔细检查, 避免出现参数错误、逻辑错误。在开始阶段, 可以利用较详尽的输出跟踪模拟运行一段时间,人工检查系统运行符合原始设定。尽可能利用模块化程序设计, 比如,在 \(\mathrm{M}/\mathrm{M}/1\) 问题模拟中,顾客到来可能遵循不同的规律,比如时齐泊松过程、非时齐泊松过程, 把产生顾客到来时刻的程序片段模块化并单独检查验证, 就可以避免在这部分出错。问题实际设定可能比较复杂, 在程序模块化以后, 如果有一种较简单的设定可以得到理论结果, 就可以用理论结果验证算法输出, 保证程序框架正确后, 再利用模块化设计修改各个模块适应实际复杂设定。

\section*{3.5 统计研究与随机模拟}

现代统计学期刊发表的论文中一半以上包括随机模拟结果 (也叫数值结果)。随机模拟可以辅助说明新模型或新方法的有效性, 尤其是在很难获得关于有效性的理论结果的情况下让我们也能说明自己的新方法的优点。

比如, 为了掌握统计量的抽样分布经常需要依靠随机模拟。在独立正态样本情况下, 我们已经有样本均值和样本方差的分布,但是对其它分布 \(F\left( x\right)\) 的样本,我们一般只能得到统计量 \(\widehat{\theta }\) 的大样本性质,在中小样本情况下很难得到 \(\widehat{\theta }\) 的抽样分布理论结果,随机模拟可以解决这样的问题。在抽取 \(F\left( x\right)\) 的很多批样本后,得到很多个统计量 \(\widehat{\theta }\) 样本值,从这些样本值可以估计 \(\widehat{\theta }\) 的抽样分布,并研究样本量多大时 \(\widehat{\theta }\) 的抽样分布与大样本极限分布相符,计算估计的均方误差, 等等。

下面举一个置信区间比较的例子说明统计研究中随机模拟的做法。实际应用中经常需要计算某个百分比的置信区间, 比如北京市成年人口中希望房价降低的人的百分比的置信区间。

设随机抽取了样本量为 \(n\) 样本,其中成功个数为 \(X\) 个,则 \(X \sim  \mathrm{B}\left( {n,p}\right) ,p\) 为成功概率,要计算 \(p\) 的置信度为 \(1 - \alpha\) 的置信区间。记样本中成功百分比为 \(\widehat{p} = X/n\) 。假设样本量 \(n \geq  {30}\) 。

当 \(n\) 较大时由中心极限定理可知

\[
W = \frac{\widehat{p} - p}{\sqrt{p\left( {1 - p}\right) }/\sqrt{n}} \tag{3.70}
\]

近似服从标准正态分布。因为 \(\widehat{p} \rightarrow  p\) a.s. \(\left( {n \rightarrow  \infty }\right)\) ,所以把(3.70)式分母中的 \(p\) 换成 \(\widehat{p}\) ,仍有

\[
Z = \frac{\widehat{p} - p}{\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right) }/\sqrt{n}}
\]

近似服从标准正态分布。于是可以构造 \(p\) 的置信度 \(1 - \alpha\) 的置信区间为

\[
\widehat{p} \pm  \frac{\lambda }{\sqrt{n}}\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right) } \tag{3.71}
\]

其中 \(\lambda  = {z}_{1 - \frac{\alpha }{2}}\) 是标准正态分布的 \(1 - \frac{\alpha }{2}\) 分位数。当 \(n\) 很大时这个置信区间应该是比较精确的,即 \(p\) 落入置信区间的实际概率应该很接近于标称的置信度 \(1 - \alpha\) ,但是在 \(n\) 不是太大的时候, \(p\) 落入置信区间的实际概率就可能与标称的置信度 \(1 - \alpha\) 有一定的差距,我们要用随机模拟考察 \(1 - \alpha ,n,p\) 的不同取值下置信区间(3.71)的覆盖率。称置信区间实际包含真实参数的概率为覆盖率。

参数 \(p\) 的另一个置信区间,称为 Wilson 置信区间,是直接求解关于 \(p\) 求解不等式

\[
\left| \frac{\widehat{p} - p}{\sqrt{p\left( {1 - p}\right) }/\sqrt{n}}\right|  \leq  \lambda  \tag{3.72}
\]

记

\[
\widetilde{p} = \frac{\widehat{p} + \frac{{\lambda }^{2}}{2n}}{1 + \frac{{\lambda }^{2}}{n}},\;\delta  = \frac{\lambda }{\sqrt{n}}\frac{\sqrt{\widehat{p}\left( {1 - \widehat{p}}\right)  + \frac{{\lambda }^{2}}{4n}}}{1 + \frac{{\lambda }^{2}}{n}},
\]

则求解(3.72)得到的 \(p\) 的置信区间为

\[
\widetilde{p} \pm  \delta  \tag{3.73}
\]

我们用随机模拟来考察(3.71)和(3.73)在不同 \(n,p\) 下的覆盖率并比较这两种置信区间。

好的置信区间应该满足如下两条要求:

(1)覆盖率大于等于置信度,两者越接近越好；

(2)置信区间越短越好。

下面设计这个模拟试验。设总共重复试验 \(M\) 次。 \(M\) 越大,对覆盖率和置信区间长度期望估计的随机误差越小。比如我们希望模拟计算的覆盖率误差控制在 0.1\% 以下 (在 95\% 置信度下)。对每组样本, 真值是否落入计算得到的置信区间是一个成败型试验, 设真实覆盖率为 \(r,M\) 次重复试验中置信区间包含真实参数的次数为 \(V\) ,则 \(V \sim  \mathrm{B}\left( {M,r}\right)\) ,覆盖率的估计值为 \(V/M\) ,此估计的标准误差为 \(\sqrt{r\left( {1 - r}\right) }/\sqrt{M}\) 。如果覆盖率 \(r\) 近似等于置信区间的置信度 \(1 - \alpha\) ,则若 \(1 - \alpha  \geq  {0.8}\) ,近似地有 \(r\left( {1 - r}\right)  \leq  {0.8} \times  {0.2} = {0.16},M\) 次试验的覆盖率估计的标准误差为 \(\sqrt{0.16}/\sqrt{M} = {0.4}/\sqrt{M}\) ,以 2 倍标准误差作为覆盖率估计的误差大小界限(根据中心极限定理可知有 \({95}\%\) 的概率使得 \(\left| {\frac{V}{M} - r}\right|\) 小于 2 倍标准误差),令 \(2 \times  {0.4}/\sqrt{M} \leq  {0.001}\) 则有 \(M = {640000}\) ,这个试验重复量很大,如果程序耗时过长我们只好降低对随机误差幅度的要求。对于更复杂的问题,如果难以得到所需重复次数 \(M\) 的理论值,可以逐步增加 \(M\) ,直到结果在需要的精度上不再变化为止。

比较 (3.71) 和 (3.73)要考虑多种 \(\left( {n,p,1 - \alpha }\right)\) 组合的影响。取 \(1 - \alpha  = {0.99},{0.95},{0.90},{0.80}\) 几种。(n, p)的值要考虑到各种不同情况,初步选取如下的一些(n, p)组合:

\[
n = {30},p = {0.1},{0.3},{0.5},{0.7},{0.9};
\]

\[
n = {120},p = {0.05},{0.1},{0.3},{0.5},{0.7},{0.9},{0.95}
\]

\[
n = {480},p = {0.01},{0.05},{0.1},{0.3},{0.5},{0.7},{0.9},{0.95},{0.99}
\]

这样,我们设计了 \(\left( {n,p,1 - \alpha }\right)\) 的 \(4 \times  \left( {5 + 7 + 9}\right)  = {84}\) 种不同情况,每种可以得到四个数值:置信区间(3.71)是否包含参数真值 (1 为包含, 0 为不包含), 置信区间(3.73)是否包含参数真值,置信区间(3.71)的长度,置信区间(3.73)的长度。对每种给定 \(\left( {n,p,1 - \alpha }\right)\) 的情况,我们重复试验 \(M\) 次,分别得到置信区间 (3.71) 的覆盖率 \({\widehat{r}}_{1}\) ,置信区间 (3.73) 的覆盖率 \({\widehat{r}}_{2}\) ,置信区间(3.71)的平均长度 \({\bar{l}}_{1}\) 和长度标准差 \({s}_{1}\) ,置信区间(3.73)的平均长度 \({\bar{l}}_{2}\) 和长度标准差 \({s}_{2}\) 。

实际编程时, 我们先编写 84 种不同情况中的一种情况去试验, 首先试验较少的重复数, 比如 \(M = {100}\) 。当程序检查无误后,再逐步增大重复次数看计算时间和存储能力是否可以接受。对此例发现 \(M = {640000}\) 可行。最后,用循环处理所有 84 种情况,每种情况重复试验 \(M = {640000}\) 次。最后的结果汇总成表3.1(限于篇幅,结果有删减)。表中每一行是一种 \(\left( {1 - \alpha ,n,p}\right)\) 的组合。注意, \({\bar{l}}_{1}\) 和 \({\bar{l}}_{2}\) 的精度可以用其标准误差 \({\mathrm{{SE}}}_{1} = {s}_{1}/\sqrt{M}\) 和 \({\mathrm{{SE}}}_{2} = {s}_{2}/\sqrt{M}\) 来估计。

从表3.1的结果看出,即使样本量已经很大 \(\left( {n = {480}}\right)\) ,公式 (3.71) 的覆盖率仍然会低于置信度 \(1 - \alpha\) ,特别是当 \(p\) 靠近 0 或 1 的情况下公式(3.71)给出的置信区间更差。另一方面,公式(3.73)的覆盖率除少数例外总是略高于标称的置信度而且与标称置信度差距一般不大,两个公式得到的置信区间的平均长度相近。但是, \(1 - \alpha  = {0.8},n = {480},p = {0.01}\) 的 Wilson 置信区间覆盖率只有 74.7\%,说明 Wilson 置信区间仍有改进必要。

\section*{3.6 Bootstrap 方法 *}

\section*{3.6.1 标准误差}

在统计建模中,伴随着参数的估计值,应该同时给出估计的 “标准误差”。设总体 \(X \sim  F\left( {x,\theta }\right) ,\theta  \in  \Theta ,\widehat{\phi }\) 是总体的一个参数 \(\phi\) 的估计量,称 \(\mathrm{{SE}} = \sqrt{\operatorname{Var}\left( \widehat{\phi }\right) }\) 为 \(\widehat{\phi }\) 的标准误差。 实际工作中 \(\mathrm{{SE}}\) 一般是未知的, \(\mathrm{{SE}}\) 的估计也称为 \(\widehat{\phi }\) 的标准误差。对有偏估计,除了标准误差外我们还希望能够估计偏差。进一步地,我们还可能希望得到统计量 \(\widehat{\phi }\) 的分布,称为抽样分布。

例 3.6.1. 设 \({X}_{i},i = 1,\ldots ,n\) 是总体 \(X \sim  F\left( x\right)\) 的样本,样本平均值 \(\widehat{\phi } = \bar{X} = \frac{1}{n}\mathop{\sum }\limits_{i}{X}_{i}\) 为 \(\phi  = {EX}\) 的点估计, \(\operatorname{SE}\left( \bar{X}\right)  = \sqrt{\operatorname{Var}\left( X\right) /n}\) ,可以用 \(S/\sqrt{n}\) 估计 \(\mathrm{{SE}}\left( {S}^{2}\right.\) 为样本方差)。 根据中心极限定理和强大数律,当样本量 \(n\) 较大时可以取 \({EX}\) 的近似 \({95}\%\) 置信区间为 \(\bar{X} \pm  2\operatorname{SE}\left( \bar{X}\right)\) 。

例 3.6.2. 考虑线性模型中参数估计的精度。设模型为

\[
\mathbf{Y} = \mathbf{X}\mathbf{\beta } + \varepsilon , \tag{3.74}
\]

其中 \(\varepsilon  \sim  \mathrm{N}\left( {0,{\sigma }^{2}{I}_{n}}\right) ,{\sigma }^{2}\) 未知, \(\mathbf{\beta }\) 是未知系数向量, \(X\) 是已知的 \(n \times  p\) 数值矩阵, \(n > p\) 。在 \(X\) 列满秩时 \(\mathbf{\beta }\) 的最小二乘估计为 \(\widehat{\mathbf{\beta }} = {\left( {X}^{T}X\right) }^{-1}{X}^{T}\mathbf{Y}\) ,而 \(\widehat{\mathbf{\beta }}\) 的协方差阵为 \(\operatorname{Var}\left( \widehat{\mathbf{\beta }}\right)  = {\sigma }^{2}{\left( {X}^{T}X\right) }^{-1}\) 。 所以,第 \(j\) 个系数 \({\beta }_{j}\) 的标准误差可估计为 \(\mathrm{{SE}}\left( {\widehat{\theta }}_{j}\right)  = \widehat{\sigma }\sqrt{{a}^{\left( jj\right) }}\) ,其中 \(\widehat{\sigma }\) 是 \(\sigma\) 的估计, \({a}^{\left( ij\right) }\) 为 \({\left( {X}^{T}X\right) }^{-1}\) 的(i, j)元素。

例 3.6.3. 设总体 \(X \sim  p\left( {x,\theta }\right) ,\theta  \in  \Theta ,{X}_{1},{X}_{2},\ldots ,{X}_{n}\) 为 \(X\) 的简单随机样本, \(\widehat{\theta }\) 是真值 \(\theta\) 的最大似然估计。在适当正则性条件下, \(\widehat{\theta }\) 渐近正态分布,渐近方差为 \(\frac{1}{n}{I}^{-1}\left( \theta \right) ,I\left( \theta \right)\) 为参数 \(\theta\) 的

表 3.1: 百分比的两种置信区间的模拟比较结果 (有删减)

\begin{center}
\adjustbox{max width=\textwidth}{
\begin{tabular}{|c|c|c|c|c|c|c|c|c|}
\hline
置信度 & 样本量 & \(p\) & \({r}_{1}\) & \({r}_{2}\) & \({\bar{l}}_{1}\) & \({\bar{l}}_{2}\) & \({s}_{1}\) & \({s}_{2}\) \\
\cline{1-9}
\multirow{16}{*}{0.99} & \multirow{3}{*}{30} & 0.1 & 95.7\% & 99.2\% & 0.2646 & 0.2870 & 0.0834 & 0.0448 \\
\cline{3-9}
 &  & 0.3 & 96.8\% & 99.2\% & 0.4221 & 0.3905 & 0.0375 & 0.0267 \\
\cline{3-9}
 &  & 0.5 & 98.4\% & 99.5\% & 0.4623 & 0.4196 & 0.0115 & 0.0084 \\
\cline{2-9}
 & \multirow{4}{*}{120} & 0.05 & 94.2\% & 99.3\% & 0.1000 & 0.1088 & 0.0204 & 0.0164 \\
\cline{3-9}
 &  & 0.1 & 98.3\% & 99.0\% & 0.1394 & 0.1423 & 0.0176 & 0.0154 \\
\cline{3-9}
 &  & 0.3 & 98.7\% & 98.8\% & 0.2144 & 0.2099 & 0.0088 & 0.0080 \\
\cline{3-9}
 &  & 0.5 & 98.7\% & 99.2\% & 0.2342 & 0.2280 & 0.0014 & 0.0013 \\
\cline{2-9}
 & \multirow{9}{*}{480} & 0.01 & 95.2\% & 99.0\% & 0.0227 & 0.0264 & 0.0057 & 0.0045 \\
\cline{3-9}
 &  & 0.05 & 98.1\% & 99.1\% & 0.0510 & 0.0521 & 0.0049 & 0.0046 \\
\cline{3-9}
 &  & 0.1 & 98.7\% & 98.8\% & 0.0703 & 0.0707 & 0.0043 & 0.0042 \\
\cline{3-9}
 &  & 0.3 & 99.0\% & 98.9\% & 0.1076 & 0.1070 & 0.0022 & 0.0021 \\
\cline{3-9}
 &  & 0.5 & 99.1\% & 99.1\% & 0.1174 & 0.1166 & 0.0002 & 0.0002 \\
\cline{3-9}
 &  & 0.7 & 99.0\% & 98.9\% & 0.1076 & 0.1070 & 0.0022 & 0.0021 \\
\cline{3-9}
 &  & 0.9 & 98.7\% & 98.8\% & 0.0703 & 0.0707 & 0.0043 & 0.0042 \\
\cline{3-9}
 &  & 0.95 & 98.1\% & 99.1\% & 0.0510 & 0.0521 & 0.0049 & 0.0046 \\
\cline{3-9}
 &  & 0.99 & 95.3\% & 99.0\% & 0.0227 & 0.0264 & 0.0057 & 0.0045 \\
\cline{1-9}
\multirow{16}{*}{0.8} & \multirow{3}{*}{30} & 0.1 & 74.3\% & 88.4\% & 0.1317 & 0.1370 & 0.0414 & 0.0323 \\
\cline{3-9}
 &  & 0.3 & 75.7\% & 84.0\% & 0.2100 & 0.2058 & 0.0186 & 0.0170 \\
\cline{3-9}
 &  & 0.5 & 80.0\% & 80.0\% & 0.2300 & 0.2241 & 0.0057 & 0.0053 \\
\cline{2-9}
 & \multirow{4}{*}{120} & 0.05 & 77.7\% & 86.4\% & 0.0498 & 0.0510 & 0.0101 & 0.0095 \\
\cline{3-9}
 &  & 0.1 & 76.8\% & 83.1\% & 0.0693 & 0.0698 & 0.0087 & 0.0084 \\
\cline{3-9}
 &  & 0.3 & 80.6\% & 80.6\% & 0.1067 & 0.1061 & 0.0044 & 0.0043 \\
\cline{3-9}
 &  & 0.5 & 76.6\% & 83.0\% & 0.1165 & 0.1157 & 0.0007 & 0.0007 \\
\cline{2-9}
 & \multirow{9}{*}{480} & 0.01 & 80.4\% & 74.7\% & 0.0113 & 0.0118 & 0.0028 & 0.0026 \\
\cline{3-9}
 &  & 0.05 & 78.8\% & 82.9\% & 0.0254 & 0.0255 & 0.0024 & 0.0024 \\
\cline{3-9}
 &  & 0.1 & 79.7\% & 80.4\% & 0.0350 & 0.0350 & 0.0021 & 0.0021 \\
\cline{3-9}
 &  & 0.3 & 80.3\% & 78.6\% & 0.0535 & 0.0535 & 0.0011 & 0.0011 \\
\cline{3-9}
 &  & 0.5 & 81.5\% & 81.5\% & 0.0584 & 0.0583 & 0.0001 & 0.0001 \\
\cline{3-9}
 &  & 0.7 & 80.4\% & 78.7\% & 0.0535 & 0.0535 & 0.0011 & 0.0011 \\
\cline{3-9}
 &  & 0.9 & 79.8\% & 80.5\% & 0.0350 & 0.0350 & 0.0021 & 0.0021 \\
\cline{3-9}
 &  & 0.95 & 78.9\% & 82.8\% & 0.0254 & 0.0255 & 0.0024 & 0.0024 \\
\cline{3-9}
 &  & 0.99 & 80.4\% & 74.6\% & 0.0113 & 0.0118 & 0.0028 & 0.0026 \\
\cline{1-9}
\hline
\end{tabular}
}
\end{center}

信息量 (参见茆诗松等 \({\left( {2006}\right) }^{\left\lbrack  7\right\rbrack  }§{2.5.2}\) 定理 2.14):

\[
I\left( \theta \right)  = \mathrm{E}\left\lbrack  {\left( \frac{\partial \ln p\left( {X,\theta }\right) }{\partial \theta }\right) }^{2}\right\rbrack   = \operatorname{Var}\left( \frac{\partial \ln p\left( {X,\theta }\right) }{\partial \theta }\right)  \tag{3.75}
\]

在加强的条件下还有

\[
I\left( \theta \right)  =  - \mathrm{E}\left( \frac{{\partial }^{2}\ln p\left( {X,\theta }\right) }{\partial {\theta }^{2}}\right)  \tag{3.76}
\]

可以用 \(\sqrt{{I}^{-1}\left( \widehat{\theta }\right) /n}\) 估计 \(\widehat{\theta }\) 的 \(\mathrm{{SE}}\) 。

例 3.6.4. 设总体 \(X \sim  p\left( {x,\mathbf{\theta }}\right) ,\mathbf{\theta } = \left( {{\theta }_{1},\ldots ,{\theta }_{m}}\right)\) ,记

\[
\mathbf{S}\left( \theta \right)  = {\nabla }_{\mathbf{\theta }}\ln p\left( {X,\mathbf{\theta }}\right)  = {\left( \frac{\partial \ln p\left( {X,\mathbf{\theta }}\right) }{\partial {\theta }_{1}},\ldots ,\frac{\partial \ln p\left( {X,\mathbf{\theta }}\right) }{\partial {\theta }_{m}}\right) }^{T}, \tag{3.77}
\]

\[
I\left( \theta \right)  = \operatorname{Var}\left( {\mathbf{S}\left( \theta \right) }\right) , \tag{3.78}
\]

称 \(I\left( \theta \right)\) 为信息量矩阵,其(i, j)元素为

\[
\operatorname{Cov}\left( {\frac{\partial \ln p\left( {X,\mathbf{\theta }}\right) }{\partial {\theta }_{i}},\frac{\partial \ln p\left( {X,\mathbf{\theta }}\right) }{\partial {\theta }_{j}}}\right)  \tag{3.79}
\]

在加强的条件下 \(I\left( \theta \right)  =  - E\left( {H\left( {X;\theta }\right) }\right) ,H\) 是 \(\ln p\left( {X,\mathbf{\theta }}\right)\) 关于自变量 \(\mathbf{\theta }\) 的海色阵,其(i, j)元素为 \(\frac{{\partial }^{2}\ln p\left( {X,\theta }\right) }{\partial {\theta }_{i}\partial {\theta }_{j}}\) 。设 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) 为 \(X\) 的简单随机样本, \(\widehat{\mathbf{\theta }}\) 为 \(\mathbf{\theta }\) 的最大似然估计,在适条件下 \(\widehat{\mathbf{\theta }}\) 渐近正态分布 \(N\left( {\mathbf{\theta },\frac{1}{n}{I}^{-1}\left( \mathbf{\theta }\right) }\right)\) ,可以用 \(- \frac{1}{n}{H}^{-1}\left( {X;\widehat{\theta }}\right)\) 作为 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 的估计。

\section*{3.6.2 Bootstrap 方法的引入}

计算参数估计的标准误差不一定总有简单的公式。例如,需要估计的参数不一定是 \({EX}\) 这样的简单特征,像中位数、相关系数这样的参数估计的标准误差就比 \({EX}\) 的估计的标准误差要困难得多。在线性模型估计的例子中, 如果独立性、线性或者正态分布的假定不满足则求参数估计方差阵变得很困难, 比如稳健回归系数的标准误差就很难得到理论公式。在最大似然估计问题中,最大似然估计不一定总是渐近正态的,信息量有时不存在或难以计算,从而无法用上面的方法给出标准误差。

设总体 \(X\) 服从某个未知分布 \(F\left( x\right) ,\mathbf{X} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 是 \(X\) 的一个样本, \(\phi\) 是 \(F\) 的一个参数,可以把 \(\phi\) 看成 \(F\) 的一个泛函 \(\phi \left( F\right)\) ,用统计量 \(\widehat{\phi } = g\left( \mathbf{X}\right)\) 估计 \(\phi\) ,设 \(\psi  = \psi \left( {g,F,n}\right)\) 是统计量 \(\widehat{\phi }\) 的某种分布特征 \((\widehat{\phi }\) 的抽样分布的数字特征)。例如 \(\psi  = \sqrt{\operatorname{Var}\left( \bar{X}\right) }\) 为统计量 \(\bar{X}\) 的标准误差,又如取 \(\psi  = E\widehat{\phi } - \phi\) 为统计量 \(\widehat{\phi }\) 的偏差。可以用随机模拟的方法估计 \(\psi\) 。

用随机模拟方法估计 \(\psi\) 的步骤如下。

( 1 )从样本 \(\mathbf{X}\) 估计总体分布 \(F\) 为 \(\widehat{F}\) ；

(2) 从 \(\widehat{F}\) 抽取 \(B\) 个独立样本 \({\mathbf{Y}}^{\left( b\right) },b = 1,\ldots ,B\) ,每一个 \({\mathbf{Y}}^{\left( b\right) }\) 样本量为 \(n\) ,称 \({\mathbf{Y}}^{\left( b\right) }\) 为 bootstrap 样本。

( 3 )从每个 bootstrap 样本 \({\mathbf{Y}}^{\left( b\right) }\) 可以估计得到 \({\widehat{\phi }}^{\left( b\right) } = g\left( {\mathbf{Y}}^{\left( b\right) }\right) ,b = 1,\ldots ,B\) 。

(4) \({\widehat{\phi }}^{\left( b\right) },b = 1,\ldots ,B\) 是 \(g\left( \mathbf{Y}\right)\) 在 \(\widehat{F}\) 下的独立同分布样本,可以用标准的估计方法估计关于 \(g\left( \mathbf{Y}\right)\) 在 \(\widehat{F}\) 下的分布特征 \(\widehat{\psi } = \psi \left( {g,\widehat{F},n}\right)\) ,估计结果记作 \(\widetilde{\psi }\) ,并以 \(\widetilde{\psi }\) 作为统计量 \(\widehat{\phi }\) 的抽样分布的数字特征 \(\psi \left( {g,F,n}\right)\) 的估计值。

从样本 \(\mathbf{X}\) 估计 \(\widehat{F}\) 时,可以采用参数模型,也可以采用经验分布函数 \({F}_{n}\) 。参数模型在模型正确时效率较高; 经验分布法使用简单,基本不依赖于模型。从经验分布 \({F}_{n}\) 抽样,相当于从 \(\mathbf{X} = \left( {{x}_{1},\ldots ,{x}_{n}}\right)\) 独立有放回抽样。

估计量的标准误差可以用 bootstrap 方法估计。

例 3.6.5. 设(H, W)为某地小学五年级学生的身高和体重的总体, \(\left( {H,W}\right)  \sim  F\left( {\cdot , \cdot  }\right)\) ,考虑 \(H\) 和 \(W\) 的相关系数 \(\phi\) 的估计。设调查了 \(n = {10}\) 个学生的身高和体重的数据 \(\left( {{h}_{i},{w}_{i}}\right) ,i = 1,2,\ldots ,n\) :

\begin{center}
\adjustbox{max width=\textwidth}{
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|c|}
\hline
\({h}_{i}\) & 144 & 166 & 163 & 143 & 152 & 169 & 130 & 159 & 160 & 175 \\
\cline{1-11}
\({w}_{i}\) & 38 & 44 & 41 & 35 & 38 & 51 & 23 & 51 & 46 & 51 \\
\cline{1-11}
\hline
\end{tabular}
}
\end{center}

计算得 \(\widehat{\phi } = g\left( {{h}_{1},{w}_{1},\ldots ,{h}_{n},{w}_{n}}\right)  = {0.904}\) 。令 \(\mathrm{{SE}}\left( \widehat{\phi }\right)  = {\left\lbrack  \operatorname{Var}\left( \widehat{\phi }\right) \right\rbrack  }^{1/2} = \psi \left( {g,F,n}\right)\) 。设 \(\widehat{F}\) 为 \(F\) 的估计,取为经验分布 \({F}_{n}\) ,则 bootstrap 方法用随机模拟方法估计 \(\psi \left( {g,{F}_{n},n}\right)\) ,然后当作 \(\operatorname{SE}\left( \widehat{\phi }\right)\) 的估计。计算步骤如下:

(1) 从 \({F}_{n}\) 中作 \(n = {10}\) 次独立抽样,即从 \(\left\{  {\left( {{h}_{1},{w}_{1}}\right) ,\ldots ,\left( {{h}_{n},{w}_{n}}\right) }\right\}\) 中有放回独立抽取 \(n\) 次, 得到 \(\widehat{F} = {F}_{n}\) 的一组样本 \({\mathbf{Y}}^{\left( 1\right) } = \left( {\left( {{h}_{1}^{\left( 1\right) },{w}_{1}^{\left( 1\right) }}\right) ,\ldots ,\left( {{h}_{n}^{\left( 1\right) },{w}_{n}^{\left( 1\right) }}\right) }\right)\) ;

(2)重复第(1)步,直到获取了 \(B\) 组 bootstrap 样本 \({\mathbf{Y}}^{\left( b\right) },b = 1,\ldots ,B\) ；

( 3 )对每一样本 \({\mathbf{Y}}^{\left( b\right) }\) 计算样本相关系数 \({\widehat{\phi }}^{\left( b\right) } = g\left( {\mathbf{Y}}^{\left( b\right) }\right) ,b = 1,\ldots ,B\) ；

(4) 把 \({\widehat{\phi }}^{\left( b\right) },b = 1,\ldots ,B\) 作为 \(\widehat{F}\) 下 \(n = {10}\) 的样本相关系数的简单随机样本,估计其样本标准差 \(S\) ,以 \(S\) 作为 \(\psi \left( {g,\widehat{F},n}\right)\) 的估计,进而用 \(S\) 估计 \(\widehat{\phi }\) 在真实的总体分布 \(F\) 下的标准误差 \(\mathrm{{SE}}\left( \widehat{\phi }\right)\) 。

取 \(B = {10000}\) 的一次 bootstrap 计算得到的标准误差估计为 \(S = {0.101}\) 。当 \(B \rightarrow  \infty\) 时 \(S \rightarrow  \psi \left( {g,{F}_{n},n}\right)\) ,但是要注意,由于抽样误差影响, \(\psi \left( {g,{F}_{n},n}\right)\) 和 \(\psi \left( {g,F,n}\right)\) 之间的误差无法避免。

也可以用参数方法估计 \(\widehat{F}\) ,比如从历史经验知道总体的身高、体重服从联合正态分布, 就可以按照联合正态总体模型从样本中得到参数最大似然估计后作为 \(\widehat{F}\) 的参数,这时 \(\widehat{F}\) 是一个参数确定的二元联合正态分布 \(\mathrm{N}\left( {{156.1},{41.8},{13.78}^{2},{8.85}^{2},{0.904}}\right)\) 。从 \(\widehat{F}\) 中独立抽样 \(n\) 个得到一组样本,共生成 \(B\) 组这样的样本,称为 bootstrap 样本。接下来的步骤只要按照上面的 (3)、(4) 估计 \(\mathrm{{SE}}\left( \widehat{\phi }\right)\) 就可以了。取 \(B = {10000}\) 的一次 bootstrap 计算得到的标准误差估计为 0.080 。

例 3.6.6. 考虑回归模型系数估计的标准误差计算。一般的回归模型可以写成

\[
{y}_{i} = h\left( {{\mathbf{x}}_{i},\mathbf{\beta }}\right)  + {\varepsilon }_{i},\;i = 1,2,\ldots ,n \tag{3.80}
\]

其中 \(h\) 已知, \(\mathbf{\beta }\) 是未知参数向量, \(\left\{  {\varepsilon }_{i}\right\}\) iid \(F\left( x\right) ,F\left( x\right)\) 未知, \(\left\{  {\mathbf{x}}_{i}\right\}\) 是确定数值向量。可以用最小二乘等方法得到 \(\mathbf{\beta }\) 的估计 \(\widehat{\mathbf{\beta }} = g\left( {{y}_{1},\ldots ,{y}_{n}}\right)\) ,希望估计参数估计的协方差阵 \(\operatorname{Var}\left( \widehat{\mathbf{\beta }}\right)\) ,协方差阵主对角线元素的平方根就是单个系数估计的标准误差。这个模型中, 未知的分布信息包括 \(\mathbf{\beta }\) 和 \(F\) 。 \(\mathbf{\beta }\) 可用 \(\widehat{\mathbf{\beta }}\) 估计, \(F\) 可以用回归残差的经验分布来估计或假设一个参数模型估计模型参数。

用 bootstrap 方法估计 \(\operatorname{Var}\left( \widehat{\mathbf{\beta }}\right)\) 的步骤如下:

( 1 )估计 \(\mathbf{\beta }\) 得到 \(\widehat{\mathbf{\beta }} = g\left( {{y}_{1},\ldots ,{y}_{n}}\right)\) ；

(2)计算残差 \({e}_{i} = {y}_{i} - h\left( {{\mathbf{x}}_{i},\widehat{\mathbf{\beta }}}\right) ,i = 1,\ldots ,n\) ；

(3) 对 \(b = 1,\ldots ,B\) 重复:从 \(\left\{  {{e}_{1},\ldots ,{e}_{n}}\right\}\) 中有放回独立抽取 \(n\) 次得 \(\left\{  {{e}_{i}^{\left( b\right) },i = 1,\ldots ,n}\right\}\) ;

(4) 令 \({y}_{i}^{\left( b\right) } = h\left( {{\mathbf{x}}_{i},\widehat{\mathbf{\beta }}}\right)  + {e}_{i}^{\left( b\right) },i = 1,\ldots ,n,b = 1,2,\ldots ,B\) ;

( 5 )对 \(b = 1,2,\ldots ,B\) 重复:从 \(\left( {{y}_{1}^{\left( b\right) },\ldots ,{y}_{n}^{\left( b\right) }}\right)\) 中估计 \({\widehat{\mathbf{\beta }}}^{\left( b\right) } = g\left( {{y}_{1}^{\left( b\right) },\ldots ,{y}_{n}^{\left( b\right) }}\right)\) 。

(6) 用 \({\widehat{\mathbf{\beta }}}^{\left( b\right) },b = 1,\ldots ,B\) 的样本方差阵估计 \(\operatorname{Var}\left( \widehat{\mathbf{\beta }}\right)\) 。

\section*{3.6.3 Bootstrap 偏差校正}

设 \(\mathbf{X} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 为总体 \(F\left( \cdot \right)\) 的样本,总体参数 \(\phi  = \phi \left( F\right)\) 的估计为 \(\widehat{\phi } = g\left( \mathbf{X}\right)\) , \(b = E\widehat{\phi } - \phi\) 为估计偏差, \(\operatorname{Var}\left( \widehat{\phi }\right)\) 为估计方差,估计的均方误差可以分解为

\[
E{\left\lbrack  \widehat{\phi } - \phi \right\rbrack  }^{2} = \operatorname{Var}\left( \widehat{\phi }\right)  + {b}^{2}. \tag{3.81}
\]

如果 \(\widehat{b}\) 是 \(b\) 的估计,则参数 \(\phi\) 的一个改善的估计为 \(\widetilde{\phi } = \widehat{\phi } - \widehat{b}\) ,新的估计在减小了偏差的同时一般也减小了均方误差。设 \(b = \psi \left( {g,F,n}\right) ,\widehat{F}\) 是总体分布 \(F\) 的一个估计,这里 \(\widehat{F}\) 取为经验分布 \({F}_{n}\) ,则可以用 \(\widehat{b} = \psi \left( {g,\widehat{F},n}\right)  = {Eg}\left( \mathbf{Y}\right)  - \widehat{\phi }\) 来估计偏差,其中 \(\mathbf{Y}\) 是总体 \({F}_{n}\) 的样本量为 \(n\) 的样本, \(\widehat{\phi }\) 恰好是总体分布为 \({F}_{n}\) 时的参数 \(\phi\) ,即 \(\widehat{\phi } = \phi \left( {F}_{n}\right)\) 。如果 \(\widehat{b}\) 不能通过理论公式计算,可以用 bootstrap 方法估计 \(\widehat{b}\) ,步骤如下:

( 1 )从 \(\left\{  {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right\}\) 独立有放回地抽取 \(n\) 个,记为 \({\mathbf{Y}}^{\left( 1\right) } = \left( {{Y}_{1}^{\left( 1\right) },{Y}_{2}^{\left( 1\right) },\ldots ,{Y}_{n}^{\left( 1\right) }}\right)\) 。

(2)重复第 (1) 步,直到获取了 \(B\) 组 bootstrap 样本 \({\mathbf{Y}}^{\left( b\right) },b = 1,2,\ldots ,B\) ；

( 3 )从每个 bootstrap 样本 \({\mathbf{Y}}^{\left( b\right) }\) 可以估计得到 \({\widehat{\phi }}^{\left( b\right) } = g\left( {\mathbf{Y}}^{\left( b\right) }\right) ,b = 1,\ldots ,B\) 。

(4) 用 \({\widehat{\phi }}^{\left( b\right) },b = 1,\ldots ,B\) 作为 \(g\left( \mathbf{Y}\right)\) 在 \({F}_{n}\) 下的独立同分布样本,估计 \(\widehat{b} = \psi \left( {g,{F}_{n},n}\right)\) 为 \(\widetilde{b} = \frac{1}{B}\mathop{\sum }\limits_{{b = 1}}^{B}{\widehat{\phi }}^{\left( b\right) } - \widehat{\phi }\) 。

最后,取 \(\widetilde{\phi } = \widehat{\phi } - \widetilde{b} = 2\widehat{\phi } - \frac{1}{B}\mathop{\sum }\limits_{{b = 1}}^{B}{\widehat{\phi }}^{\left( b\right) }\) 为改善的估计。

例 3.6.7. 设 \(X \sim  \mathrm{N}\left( {\mu ,{\sigma }^{2}}\right) ,\mu ,{\sigma }^{2}\) 未知, \({x}_{1},{x}_{2},\ldots ,{x}_{n}\) 为 \(X\) 的样本。考虑 \(\phi  = {\mu }^{2}\) 的估计。 用最大似然估计法估计 \(\phi\) 为 \(\widehat{\phi } = {\bar{X}}^{2}\) ,其中 \(\bar{X}\) 为样本平均值。令 \(Z = \sqrt{n}\left( {\bar{X} - \mu }\right) /\sigma\) ,则 \(Z \sim  \mathrm{N}\left( {0,1}\right)\) 。可以计算出估计偏差为

\[
b = E{\bar{X}}^{2} - {\mu }^{2} = E{\left( \mu  + \frac{\sigma }{\sqrt{n}}Z\right) }^{2} - {\mu }^{2} = \frac{{\sigma }^{2}}{n}, \tag{3.82}
\]

估计的均方误差为

\[
{L}_{0} = E{\left( {\bar{X}}^{2} - {\mu }^{2}\right) }^{2} = E{\left( \frac{2\sigma \mu }{\sqrt{n}}Z + \frac{{\sigma }^{2}}{n}{Z}^{2}\right) }^{2}
\]

\[
= \frac{4{\sigma }^{2}{\mu }^{2}}{n} + \frac{3{\sigma }^{4}}{{n}^{2}}. \tag{3.83}
\]

估计 \(b\) 为 \({\widehat{b}}_{1} = {S}^{2}/n\left( {S}^{2}\right.\) 为样本方差),用 \({\widehat{\phi }}_{1} = {\bar{X}}^{2} - {\widehat{b}}_{1} = {\bar{X}}^{2} - \frac{{S}^{2}}{n}\) 作为 \({\mu }^{2}\) 的改善的估计,则 \({\widehat{\phi }}_{1}\) 的均方误差比 \(\widehat{\phi }\) 的均方误差减小了 \(\frac{n - 3}{n - 1}\frac{{\sigma }^{4}}{{n}^{2}}\) (设 \(n > 3\) ,见习题17)。

如果模型更为复杂,比如,总体分布类型未知, \({\widehat{b}}_{1}\) 这样的简单偏差估计很难得到,这种情况下可以用 bootstrap 方法进行偏差校正, 步骤如下:

(1) 对 \(b = 1,2,\ldots ,B\) 重复: 从 \({x}_{1},{x}_{2},\ldots ,{x}_{n}\) 独立有放回地抽取 \(n\) 个,组成 bootstrap 样本 \({\mathbf{Y}}^{\left( b\right) } = \left( {{y}_{1}^{\left( b\right) },\ldots ,{y}_{n}^{\left( b\right) }}\right) ;\)

(2)对每个 bootstrap 样本计算 \({\widehat{\phi }}^{\left( b\right) } = {\left( \frac{1}{n}\mathop{\sum }\limits_{{i = 1}}^{n}{y}_{i}^{\left( b\right) }\right) }^{2}\) ；

( 3 )用 \(\widetilde{\phi } = 2{\bar{X}}^{2} - \frac{1}{B}\mathop{\sum }\limits_{{b = 1}}^{B}{\widehat{\phi }}^{\left( b\right) }\) 作为 \({\mu }^{2}\) 的改善的估计。

Jacknife 方法是另外一种对估计量的偏差和方差进行估计的方法, 这种方法不需要从原来的样本重新随机抽样,而是把原来的 \(n\) 个样本点分为 \(r\) 份,每次删去其中一份后计算统计量值,利用 \(r\) 个这样的统计量值对估计量的偏差和方差进行估计。详见 \(\operatorname{Gentle}{\left( {2002}\right) }^{\left\lbrack  {18}\right\rbrack  }§{3.3}\) 。

\section*{3.6.4 Bootstrap 置信区间}

枢轴量法是构造置信区间的最基本的方法。设 \(\phi\) 是总体 \(F\left( \cdot \right)\) 的一个参数,看成 \(F\) 的一个泛函 \(\phi  = \phi \left( F\right)\) 。 \(\mathbf{X} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 为总体的样本, \(g\left( \mathbf{X}\right)\) 为与 \(\phi\) 有关系的一个统计量, 经常是 \(\phi\) 的估计量。如果有变换 \(W = h\left( {g\left( \mathbf{X}\right) ,\phi }\right)\) 使得 \(W\) 的分布不依赖于任何未知参数,则设 \(W\) 的左右两侧分位数分别为 \({w}_{\frac{\alpha }{2}}\) 和 \({w}_{1 - \frac{\alpha }{2}}\) ,有

\[
P\left( {{w}_{\frac{\alpha }{2}} < h\left( {T,\phi }\right)  < {w}_{1 - \frac{\alpha }{2}}}\right)  = 1 - \alpha , \tag{3.84}
\]

反解上面的不等式可以得到 \(\phi\) 的置信区间。

如果对枢轴量 \(W\) 很难求分位数时,可以用 bootstrap 方法获得置信区间。设 \(\widehat{F}\) 为总体分布 \(F\) 的估计,设 \(\mathbf{Y} = \left( {{y}_{1},\ldots ,{y}_{n}}\right)\) 为总体 \(\widehat{F}\) 的样本, \(\widehat{\phi } = \phi \left( \widehat{F}\right)\) 为与总体 \(\widehat{F}\) 对应的参数 \(\phi\) 的值,实际是 \(\phi\) 的估计值,则 \(V = h\left( {g\left( \mathbf{Y}\right) ,\widehat{\phi }}\right)\) 与 \(W\) 的分布相近,可以用 \(V\) 的分位数近似 \(W\) 的分位数。

例 3.6.8. 设总体 \(X \sim  F\left( {x,\theta }\right) ,\theta\) 为总体的未知参数, \(\phi  = {EX},\mathbf{X} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 为总体的样本,则 \(g\left( \mathbf{X}\right)  = \bar{X}\) 是 \(\phi\) 的估计,若 \(W = h\left( {g\left( \mathbf{X}\right) ,\phi }\right)  = \bar{X} - {EX}\) 的分布与 \(\theta\) 无关,求 \(W\) 的分位数 \({w}_{\frac{\alpha }{2}}\) 和 \({w}_{1 - \frac{\alpha }{2}}\) 就可以构造 \(\phi  = {EX}\) 的置信区间 \(\left( {\bar{X} - {w}_{1 - \frac{\alpha }{2}},\bar{X} - {w}_{\frac{\alpha }{2}}}\right)\) 。

若 \(W\) 的分位数无法求得,用经验分布 \({F}_{n}\) 作为总体分布 \(F\) 的估计,这时 \(\phi \left( {F}_{n}\right)  = \bar{X}\) ,设 \(\mathbf{Y} = \left( {{Y}_{1},\ldots ,{Y}_{n}}\right)\) 为 \({F}_{n}\) 的样本, \(V = h\left( {g\left( \mathbf{Y}\right) ,\bar{X}}\right)  = \bar{Y} - \bar{X}\) ,这里 \(\bar{X}\) 作为已知值,可以用 \(V\) 的分位数近似代替 \(W\) 的分位数。求 \(V\) 的分位数,只要用有放回独立抽样方法从 \({x}_{1},{x}_{2},\ldots ,{x}_{n}\) 抽取 \({F}_{n}\) 的 \(B\) 组样本 \({\mathbf{Y}}^{\left( b\right) } = \left( {{y}_{1}^{\left( b\right) },\ldots ,{y}_{n}^{\left( b\right) }}\right) ,b = 1,2,\ldots ,B\) ,对每组样本计算平均值 \({\bar{Y}}^{\left( b\right) }\) ,定义 \({V}^{\left( b\right) } = {\bar{Y}}^{\left( b\right) } - \bar{X}\) ,用 \({V}^{\left( b\right) },b = 1,2,\ldots ,B\) 的样本分位数估计 \(V\) 的分位数,作为 \(W\) 的分位数 \({w}_{\frac{\alpha }{2}}\) 和 \({w}_{1 - \frac{\alpha }{2}}\) 的近似。

\section*{3.7 MCMC}

\section*{3.7.1 马氏链和 MCMC 介绍}

实际工作中经常遇到分布复杂的高维随机向量抽样问题。§3.2.4的重要抽样法可以应付维数不太高的情况, 但是对于维数很高而且分布很复杂 (比如, 分布密度多峰而且位置不易确定的情况)则难以处理。

MCMC(马氏链蒙特卡洛) 是一种对高维随机向量抽样的方法, 此方法模拟一个马氏链, 使马氏链的平稳分布为目标分布, 由此产生大量的近似服从目标分布的样本, 但样本不是相互独立的。MCMC 的目标分布密度函数或概率函数可以只计算到差一个常数倍的值。MCMC 方法适用范围广, 近年来获得了广泛的应用。

先介绍马氏链的概念。

设 \(\left\{  {{X}_{t},t = 0,1,\ldots }\right\}\) 为随机变量序列,称为一个随机过程。称 \({X}_{t}\) 为 “系统在时刻 \(t\) 的状态”。为讨论简单起见,设所有 \({X}_{t}\) 均取值于有限集合 \(S = \{ 1,2,\ldots ,m\}\) ,称 \(S\) 为状态空间。 如果 \(\left\{  {X}_{t}\right\}\) 满足

\[
P\left( {{X}_{t + 1} = j \mid  {X}_{0} = {k}_{0},\ldots {X}_{t - 1} = {k}_{t - 1},{X}_{t} = i}\right)
\]

\[
= P\left( {{X}_{t + 1} = j \mid  {X}_{t} = i}\right)  = {p}_{ij},t = 0,1,\ldots ,{k}_{0},\ldots ,{k}_{t - 1},i,j \in  S, \tag{3.85}
\]

则称 \(\left\{  {X}_{t}\right\}\) 为马氏链, \({p}_{ij}\) 为转移概率,矩阵 \(P = {\left( {p}_{ij}\right) }_{m \times  m}\) 为转移概率矩阵。显然 \(\mathop{\sum }\limits_{{j = 1}}^{m}{p}_{ij} =\)  \(1,i = 1,2,\ldots ,m\) 。对马氏链, \(P\left( {{X}_{t + k} = j \mid  {X}_{t} = i}\right) \overset{\bigtriangleup }{ = }{p}_{ij}^{\left( k\right) }\) 也不依赖于 \(t\) ,称为 \(k\) 步转移概率。 如果对任意 \(i,j \in  S,i \neq  j\) 都存在 \(k \geq  1\) 使得 \({p}_{ij}^{\left( k\right) } > 0\) 则称 \(\left\{  {X}_{t}\right\}\) 为不可约马氏链。不可约马氏链的所有状态是互相连通的,即总能经过若干步后互相转移。对马氏链 \(\left\{  {X}_{t}\right\}\) 的某个状态 \(i\) ,如果存在 \(k \geq  0\) 使得 \({p}_{ii}^{\left( k\right) } > 0\) 并且 \({p}_{ii}^{\left( k + 1\right) } > 0\) ,则称 \(i\) 是非周期的。如果一个马氏链所有状态都是非周期的, 则该马氏链称为非周期的。不可约马氏链只要有一个状态是非周期的则所有状态是非周期的。对只有有限个状态的非周期不可约马氏链有

\[
\mathop{\lim }\limits_{{n \rightarrow  \infty }}P\left( {{X}_{n} = j}\right)  = {\pi }_{j},j = 1,2,\ldots ,m, \tag{3.86}
\]

其中 \(\left\{  {{\pi }_{j},j = 1,2,\ldots ,m}\right\}\) 为常数,称为 \(\left\{  {X}_{t}\right\}\) 的极限分布。 \(\left\{  {\pi }_{j}\right\}\) 满足方程组

\[
\left\{  \begin{array}{l} \mathop{\sum }\limits_{{i = 1}}^{m}{\pi }_{i}{p}_{ij} = {\pi }_{j},j = 1,2,\ldots ,m \\  \mathop{\sum }\limits_{{j = 1}}^{m}{\pi }_{j} = 1, \end{array}\right.  \tag{3.87}
\]

称满足(3.87)的分布 \(\left\{  {\pi }_{j}\right\}\) 为平稳分布。对只有有限个状态的非周期不可约马氏链,极限分布和平稳分布存在且为同一分布。

如果允许状态空间 \(S\) 为可列个元素,极限分布的条件需要更多的讨论。对状态 \(i\) ,如从状态 \(i\) 出发总能再返回状态 \(i\) ,则称状态 \(i\) 是常返的 (recurrent)。对常返状态 \(i\) ,如果从 \(i\) 出发首次返回 \(i\) 的时间的期望有限,称 \(i\) 是正常返的。非周期正常返状态称为遍历的。所有状态都遍历的马氏链称为遍历马氏链。非周期遍历马氏链存在唯一的极限分布和平稳分布, 且二者相同。

如果存在 \(\left\{  {{\pi }_{j},j = 1,2,\ldots ,m}\right\}  ,{\pi }_{j} \geq  0,\mathop{\sum }\limits_{{j = 1}}^{m}{\pi }_{j} = 1\) ,使得

\[
{\pi }_{i}{p}_{ij} = {\pi }_{j}{p}_{ji},\forall i \neq  j, \tag{3.88}
\]

称这样的马氏链为细致平衡的 (detailed balanced),这时 \(\left\{  {\pi }_{j}\right\}\) 是 \(\left\{  {X}_{t}\right\}\) 的平稳分布。事实上, 若 \(P\left( {{X}_{t} = i}\right)  = {\pi }_{i},i \in  S\) ,则

\[
P\left( {{X}_{t + 1} = j}\right)  = \mathop{\sum }\limits_{i}{\pi }_{i}{p}_{ij} = \mathop{\sum }\limits_{i}{\pi }_{j}{p}_{ji} = {\pi }_{j}\mathop{\sum }\limits_{i}{p}_{ji} = {\pi }_{j},\forall j \in  S. \tag{3.89}
\]

马氏链的概念可以推广到 \({X}_{t}\) 的取值集合 \(\mathcal{X}\) 为可列集或 \({\mathbb{R}}^{d}\) 的区域的情形。如果各 \(\left\{  {X}_{t}\right\}\) 的有限维分布是连续型的,则 (3.85) 可以改用条件密度表示,这时的 \(\left\{  {X}_{t}\right\}\) 按照随机过程论中的习惯应该称作马氏过程, 但这里还是叫做马氏链。

如果遍历的不可约马氏链 \(\left\{  {X}_{t}\right\}\) 有平稳分布 \(\pi \left( x\right) ,x \in  \mathcal{X}\) ,则从任意初值出发模拟产生序列 \(\left\{  {X}_{t}\right\}\) ,当 \(t\) 很大时, \({X}_{t}\) 的分布就近似服从 \(\pi\) ,抛弃开始的一段后的 \({X}_{t}\) 序列可以作为分布 \(\pi\) 的相关的样本,抛弃的一段序列叫做老化期。设 \(Y \sim  \pi \left( \cdot \right) ,h\left( y\right) ,y \in  \mathcal{X}\) 是有界函数,为估计 \(\theta  \triangleq  {Eh}\left( Y\right)\) ,用 \({X}_{t},t = k + 1,\ldots ,n\) 作为 \(\pi\) 的样本,用估计量 \(\widehat{\theta } = \frac{1}{n - k}\mathop{\sum }\limits_{{t = k + 1}}^{n}h\left( {X}_{t}\right)\) 来估计 \(\theta\) ,则 \(\widehat{\theta }\) 是 \(\theta\) 的强相合估计。老化期长度 \(k\) 可以从 \(\widehat{\theta }\) 的变化图形经验地选取。

这样的估计量 \(\widehat{\theta }\) 是相关样本的平均值,无法用原来独立样本的公式估计 \(\operatorname{Var}\left( \widehat{\theta }\right)\) 从而得到 \(\widehat{\theta }\) 的标准误差。为了估计 \(\operatorname{Var}\left( \widehat{\theta }\right)\) ,可以采用如下的分段平均法。把样本 \({X}_{k + 1},\ldots ,{X}_{n}\) 分为 \(s\) 段,每段 \(r\) 个 (设 \(n - k = {sr}\) )。设第 \(j\) 段的 \(r\) 个 \(h\left( {X}_{t}\right)\) 的平均值为 \({Z}_{j},j = 1,2,\ldots ,s\) ,设 \(\left\{  {{Z}_{j},j = 1,2,\ldots ,s}\right\}\) 的样本方差为 \({\widehat{\sigma }}^{2} = \frac{1}{s - 1}\mathop{\sum }\limits_{{j = 1}}^{s}{\left( {Z}_{j} - \bar{Z}\right) }^{2}\) ,因为 \(\widehat{\theta }\) 等于 \(\left\{  {{Z}_{j},j = 1,2,\ldots ,s}\right\}\) 的样本均值 \(\bar{Z}\) ,当 \(r\) 足够大时,可以认为各 \({Z}_{j}\) 相关性已经很弱,这时 \(\widehat{\theta }\) 的方差可以用 \({\widehat{\sigma }}^{2}/s\) 估计。 \(r\) 的大小依赖于不同时刻的 \({X}_{t}\) 的相关性强弱,相关性越强,需要的 \(r\) 越大。

以上的方法就是 MCMC 方法 (马氏链蒙特卡洛)。一般地,对高维或取值空间 \(\mathcal{X}\) 结构复杂的随机向量 \(X\) , MCMC 方法构造取值于 \(\mathcal{X}\) 的马氏链,使其平稳分布为 \(X\) 的目标分布。 模拟此马氏链,抛弃开始的部分抽样值,把剩余部分作为 \(X\) 的非独立抽样。非独立抽样的估计效率比独立抽样低。

MCMC 方法的关键在于如何从第 \(t\) 时刻转移到第 \(t + 1\) 时刻。好的转移算法应该使得马氏链比较快地收敛到平稳分布,并且不会长时间地停留在取值空间 \(\mathcal{X}\) 的局部区域内(在目标分布是多峰分布且峰高度差异较大时容易出现这种问题)。

Metropolis-Hasting 方法 (MH 方法) 是一个基本的 MCMC 算法, 此算法在每一步试探地进行转移 (如随机游动),如果转移后能提高状态 \({x}_{t}\) 在目标分布 \(\pi\) 中的密度值则接受转移结果, 否则以一定的概率决定是转移还是停留不动。

Gibbs 抽样是另外一种常用的 MCMC 方法, 此方法轮流延各个坐标轴方向转移, 且转移概率由当前状态下用其它坐标预测转移方向坐标的条件分布给出。因为利用了目标分布的条件分布,所以 Gibbs 抽样方法的效率比 MH 方法效率更高。

\section*{3.7.2 Metropolis-Hasting 抽样}

设随机变量 \(X\) 分布为 \(\pi \left( x\right) ,x \in  \mathcal{X}\) 。为论述简单起见仍假设 \(\mathcal{X}\) 是离散集合。算法需要一个试转移概率函数 \(T\left( {y \mid  x}\right) ,x,y \in  \mathcal{X}\) ,满足 \(0 \leq  T\left( {y \mid  x}\right)  \leq  1,\mathop{\sum }\limits_{y}T\left( {y \mid  x}\right)  = 1\) ,并且

\[
T\left( {y \mid  x}\right)  > 0 \Leftrightarrow  T\left( {x \mid  y}\right)  > 0. \tag{3.90}
\]

算法首先从 \(\mathcal{X}\) 中任意取初值 \({X}^{\left( 0\right) }\) 。设经过 \(t\) 步后算法的当前状态为 \({X}^{\left( t\right) }\) ,则下一步由试转移分布 \(T\left( {y \mid  {X}^{\left( t\right) }}\right)\) 抽取 \(Y\) ,并生成 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,然后按如下规则转移:

\[
{X}^{\left( t + 1\right) } = \left\{  \begin{array}{ll} Y & \text{ 若 }U \leq  r\left( {{X}^{\left( t\right) },Y}\right) \\  {X}^{\left( t\right) } & \text{ 否则 } \end{array}\right.  \tag{3.91}
\]

其中

\[
r\left( {x,y}\right)  = \min \left\{  {1,\frac{\pi \left( y\right) T\left( {x \mid  y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }}\right\}  . \tag{3.92}
\]

在 MH 算法中如果取 \(T\left( {y \mid  x}\right)  = T\left( {x \mid  y}\right)\) ,则 \(r\left( {x,y}\right)  = \min \left( {1,\frac{\pi \left( y\right) }{\pi \left( x\right) }}\right)\) ,相应的算法称为 Metropolis 抽样法。

如果取 \(T\left( {y \mid  x}\right)  = g\left( y\right)\) (不依赖于 \(x\) ),则 \(r\left( {x,y}\right)  = \min \left( {1,\frac{\pi \left( y\right) /g\left( y\right) }{\pi \left( x\right) /g\left( x\right) }}\right)\) ,相应的算法称为 Metropolis 独立抽样法,和重要抽样有相似之处,试抽样分布 \(g\left( \cdot \right)\) 经常取为相对重尾的分布。

在 MH 算法中,目标分布 \(\pi \left( x\right)\) 可以用差一个常数倍的 \(\widetilde{\pi }\left( x\right)  = {C\pi }\left( x\right)\) 代替,这样关于目标分布仅知道差一个常数倍的 \(\widetilde{\pi }\left( x\right)\) 的情形,也可以使用此算法。

下面说明 \(\mathrm{{MH}}\) 抽样方法的合理性。我们来验证 \(\mathrm{{MH}}\) 抽样的转移概率 \(A\left( {x,y}\right)  = P\left( {{X}^{\left( t + 1\right) } = }\right.\)  \(\left. {y \mid  {X}^{\left( t\right) } = x}\right)\) 满足细致平衡条件。易见

\[
A\left( {x,y}\right)  = \left\{  \begin{array}{ll} T\left( {y \mid  x}\right) r\left( {x,y}\right) , & y \neq  x, \\  T\left( {x \mid  x}\right)  + \mathop{\sum }\limits_{{z \neq  x}}T\left( {z \mid  x}\right) \left\lbrack  {1 - r\left( {x,z}\right) }\right\rbrack  , & y = x, \end{array}\right.  \tag{3.93}
\]

于是当 \(x \neq  y\) 时

\[
\pi \left( x\right) A\left( {x,y}\right)  = \pi \left( x\right) T\left( {y \mid  x}\right) \min \left\{  {1,\frac{\pi \left( y\right) T\left( {x \mid  y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }}\right\}   = \min \{ \pi \left( x\right) T\left( {y \mid  x}\right) ,\pi \left( y\right) T\left( {x \mid  y}\right) \} ,
\]

等式右侧关于 \(x,y\) 是对称的,所以等式左侧把 \(x,y\) 交换后仍相等。所以, \(\mathrm{{MH}}\) 构造的马氏链以 \(\{ \pi \left( x\right) \}\) 为平稳分布。多数情况下 MH 构造的马氏链也以 \(\{ \pi \left( x\right) \}\) 为极限分布。

(3.92) 中的 \(r\left( {x,y}\right)\) 还可以推广为如下的形式

\[
\widetilde{r}\left( {x,y}\right)  = \frac{\alpha \left( {x,y}\right) }{\pi \left( x\right) T\left( {y \mid  x}\right) }, \tag{3.94}
\]

其中 \(\alpha \left( {x,y}\right)\) 是任意的满足 \(\alpha \left( {x,y}\right)  = \alpha \left( {y,x}\right)\) 且使得 \(\widetilde{r}\left( {x,y}\right)  \leq  1\) 的函数。易见这样的 \(\widetilde{r}\left( {x,y}\right)\) 仍使得生成的马氏链满足细致平衡条件。

例 3.7.1. \(X\) 的取值集合 \(\mathcal{X}\) 可能是很大的,以至于无法穷举,目标分布 \(\pi \left( x\right)\) 可能是只能确定到差一个常数倍。

例如, 设

\[
\mathcal{X} = \left\{  {\mathbf{x} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)  : \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right) \text{ 为 }\left( {1,2,\ldots ,n}\right) \text{ 的一个排列,并满足 }\mathop{\sum }\limits_{{j = 1}}^{n}j{x}_{j} > a}\right\}  ,
\]

其中 \(a\) 是一个给定的常数。用 \(\left| \mathcal{X}\right|\) 表示 \(\mathcal{X}\) 的元素个数,当 \(n\) 较大时 \(\mathcal{X}\) 是 \(\left( {1,2,\ldots ,n}\right)\) 的所有 \(n!\) 个排列的一个子集, \(\left| \mathcal{X}\right|\) 很大,很难穷举 \(\mathcal{X}\) 的元素,从而 \(\left| \mathcal{X}\right|\) 未知。

设 \(X\) 服从 \(\mathcal{X}\) 上的均匀分布,即 \(\pi \left( \mathbf{x}\right)  = C,x \in  \mathcal{X},C = 1/\left| \mathcal{X}\right|\) 但 \(C\) 未知。要用 MH 方法产生 \(X\) 的抽样序列。

试抽样 \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)\) 如果允许转移到所有的 \(\mathbf{y}\) 是很难执行的,因为 \(\mathbf{y}\) 的个数太多了。我们定义一个 \(\mathbf{x}\) 的近邻的概念,仅考虑转移到 \(\mathbf{x}\) 的近邻。一种定义是,如果把 \(\mathbf{x}\) 的 \(n\) 个元素中的某两个交换位置后可以得到 \(\mathbf{y} \in  \mathcal{X}\) ,则 \(\mathbf{y}\) 称为 \(\mathbf{x}\) 的一个近邻,记 \(N\left( \mathbf{x}\right)\) 为 \(\mathbf{x}\) 的所有近邻的集合,记 \(\left| {N\left( \mathbf{x}\right) }\right|\) 为 \(\mathbf{x}\) 的近邻的个数。当 \(n\) 很大时,求 \(N\left( \mathbf{x}\right)\) 也需要从 \({C}_{n}^{2} = \frac{1}{2}n\left( {n - 1}\right)\) 个可能的元素中用穷举法选择。取试转移概率函数为

\[
T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = \frac{1}{\left| N\left( \mathbf{x}\right) \right| },\mathbf{x},\mathbf{y} \in  \mathcal{X},
\]

即从 \(x\) 出发,等可能地试转移到 \(x\) 的任何一个近邻上。

因为目标分布 \(\pi \left( \mathbf{x}\right)\) 是常数,所以这时

\[
r\left( {\mathbf{x},\mathbf{y}}\right)  = \min \left( {1,\frac{\left| N\left( \mathbf{x}\right) \right| }{\left| N\left( \mathbf{y}\right) \right| }}\right) ,
\]

即从 \(\mathbf{x}\) 试转移到 \(\mathbf{y}\) 后,如果 \(\mathbf{y}\) 的近邻数不超过 \(\mathbf{x}\) 的近邻数则确定转移到 \(\mathbf{y}\) ,否则,仅按概率 \(\left| {N\left( \mathbf{x}\right) }\right| /\left| {N\left( \mathbf{y}\right) }\right|\) 转移到 \(\mathbf{y}\) 。这就构成了对 \(X\) 抽样的 \(\mathrm{{MH}}\) 算法。 连续型分布的 \(\mathbf{{MH}}\) 抽样法 对于连续型的目标分布,设 \(\pi \left( x\right)\) 为目标分布的密度,这时 \(T\left( {y \mid  x}\right)\) 改为给定 \(x\) 条件下的试抽样密度, \(r\left( {x,y}\right)\) 定义不变,算法和离散型目标分布的情形相同。

例 3.7.2. 考虑一个贝叶斯推断问题。在金融投资中, 投资者经常把若干种证券组合在一起来减少风险。假设有 5 支股票的 \(n = {250}\) 个交易日的收益率记录,每个交易日都找出这 5 支股票收益率最高的一个,设 \({X}_{i}\) 表示第 \(i\) 支股票在 \(n\) 个交易日中收益率为最高的次数 \(\left( {i = 1,2,\ldots ,5}\right)\) 。设 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 服从多项分布,相应的概率假设为

\[
\mathbf{p} = \left( {\frac{1}{3},\frac{1 - \beta }{3},\frac{1 - {2\beta }}{3},\frac{2\beta }{3},\frac{\beta }{3}}\right) ,
\]

其中 \(\beta  \in  \left( {0,{0.5}}\right)\) 为未知参数。假设 \(\beta\) 有先验分布 \({p}_{0}\left( \beta \right)  \sim  \mathrm{U}\left( {0,{0.5}}\right)\) 。设 \(\left( {{x}_{1},\ldots ,{x}_{5}}\right)\) 为 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 的观测值,则 \(\beta\) 的后验分布为

\[
f\left( {\beta  \mid  {x}_{1},\ldots ,{x}_{5}}\right)  \propto  p\left( {{x}_{1},\ldots ,{x}_{5} \mid  \beta }\right) {p}_{0}\left( \beta \right)
\]

\[
= \left( \begin{matrix} n \\  {x}_{1},\ldots ,{x}_{5} \end{matrix}\right) {\left( \frac{1}{3}\right) }^{{x}_{1}}{\left( \frac{1 - \beta }{3}\right) }^{{x}_{2}}{\left( \frac{1 - {2\beta }}{3}\right) }^{{x}_{3}}{\left( \frac{2\beta }{3}\right) }^{{x}_{4}}{\left( \frac{\beta }{3}\right) }^{{x}_{5}}\frac{1}{0.5}{I}_{\left( 0,{0.5}\right) }\left( \beta \right)
\]

\[
\propto  {\left( 1 - \beta \right) }^{{x}_{2}}{\left( 1 - 2\beta \right) }^{{x}_{3}}{\beta }^{{x}_{4} + {x}_{5}}{I}_{\left( 0,{0.5}\right) }\left( \beta \right) \overset{ \bigtriangleup  }{ = }\widetilde{\pi }\left( \beta \right) .
\]

为了求 \(\beta\) 后验均值,需要产生服从 \(f\left( {\beta  \mid  {x}_{1},\ldots ,{x}_{5}}\right)\) 的抽样。从 \(\beta\) 的后验分布很难直接抽样,采用 Metropolis 抽样法。设当前 \(\beta\) 的状态为 \({\beta }^{\left( t\right) }\) ,取试抽样分布 \(T\left( {y \mid  {\beta }^{\left( t\right) }}\right)\) 为 \(\mathrm{U}\left( {0,{0.5}}\right)\) , 则 \(T\left( {y \mid  x}\right)  = T\left( {x \mid  y}\right)\) ,

\[
r\left( {{\beta }^{\left( t\right) },y}\right)  = \min \left( {1,\frac{\widetilde{\pi }\left( y\right) }{\widetilde{\pi }\left( {\beta }^{\left( t\right) }\right) }}\right)  = \min \left( {1,{\left( \frac{1 - y}{1 - {\beta }^{\left( t\right) }}\right) }^{{x}_{2}}{\left( \frac{1 - {2y}}{1 - 2{\beta }^{\left( t\right) }}\right) }^{{x}_{3}}{\left( \frac{y}{{\beta }^{\left( t\right) }}\right) }^{{x}_{4} + {x}_{5}}}\right) ,
\]

从 \(\mathrm{U}\left( {0,{0.5}}\right)\) 试抽取 \(y\) ,以概率 \(r\left( {{\beta }^{\left( t\right) },y}\right)\) 接受 \({\beta }^{\left( t + 1\right) } = y\) 即可。

随机游动 \(\mathbf{{MH}}\) 算法 \({MH}\) 抽样中试转移概率函数 \(T\left( {y \mid  x}\right)\) 较难找到,容易想到的是从 \({x}^{\left( t\right) }\) 作随机游动的试转移方法, 叫做随机游动 Metropolis 抽样。

设 \(X\) 的目标分布 \(\pi \left( \mathbf{x}\right)\) 取值于欧式空间 \(\mathcal{X} = {\mathbb{R}}^{d}\) 。从 \({\mathbf{x}}^{\left( t\right) }\) 出发试转移,令

\[
\mathbf{y} = {\mathbf{x}}^{\left( t\right) } + {\mathbf{\varepsilon }}_{t}, \tag{3.95}
\]

其中 \({\varepsilon }_{t} \sim  g\left( {\mathbf{x};\sigma }\right)\) 对不同 \(t\) 是独立同分布的, \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = g\left( {\mathbf{y} - \mathbf{x}}\right)\) 。设 \(g\) 是关于 \(\mathbf{x} = \mathbf{0}\) 对称的分布,则 \(T\left( {\mathbf{y} \mid  \mathbf{x}}\right)  = T\left( {\mathbf{x} \mid  \mathbf{y}}\right)\) 。

常取 \(g\) 为 \(\mathrm{N}\left( {\mathbf{0},{\sigma }^{2}I}\right)\) 和半径为 \(\sigma\) 的中心为0的球内的均匀分布。

转移法则为: 从 \({\mathbf{x}}^{\left( t\right) }\) 出发试转移到 \(\mathbf{y}\) 后,若 \(\pi \left( \mathbf{y}\right)  > \pi \left( {\mathbf{x}}^{\left( t\right) }\right)\) 则令 \({\mathbf{x}}^{\left( t + 1\right) } = \mathbf{y}\) ; 否则,独立地抽取 \(U \sim  \mathrm{U}\left( {0,1}\right)\) ,取

\[
{\mathbf{x}}^{\left( t + 1\right) } = \left\{  \begin{array}{ll} \mathbf{y}, & \text{ 当 }U \leq  \pi \left( \mathbf{y}\right) /\pi \left( {\mathbf{x}}^{\left( t\right) }\right) , \\  {\mathbf{x}}^{\left( t\right) }, & \text{ 否则. } \end{array}\right.
\]

随机游动 \(\mathrm{{MH}}\) 算法是一种 Metropolis 抽样方法。随机游动的步幅 \(\sigma\) 是重要参数,步幅过大导致拒绝率大, 步幅过小使得序列的相关性太强, 收敛到平衡态速度太慢。一个建议选法是试验各种选法, 使得试抽样被接受的概率在 0.25 到 0.35 之间。

例 3.7.3. 考虑如下的简单气体模型: 在平面区域 \(G = \left\lbrack  {0,A}\right\rbrack   \times  \left\lbrack  {0,B}\right\rbrack\) 内有 \(K\) 个直径为 \(d\) 的刚性圆盘。随机向量 \(\mathbf{X} = \left( {{x}_{1},{y}_{1},\ldots ,{x}_{K},{y}_{k}}\right)\) 为这些圆盘的位置坐标。分布 \(\pi \left( \mathbf{x}\right)\) 是 \(G\) 内所有允许位置的均匀分布。希望对 \(\pi\) 抽样。

先找一个初始的允许位置 \({\mathbf{x}}^{\left( 0\right) }\) 。比如,把圆盘整齐地排列在左上角。

设已得到 \({\mathbf{x}}^{\left( t\right) }\) ,随机选取一个圆盘 \(i\) ,把圆盘 \(i\) 的位置试移动到 \(\left( {{x}_{i}^{\prime },{y}_{i}^{\prime }}\right)  = \left( {{x}_{i} + {\delta }_{i},{y}_{i} + {\epsilon }_{i}}\right)\) , 其中 \({\delta }_{i},{\epsilon }_{i}\) 独立同 \(\mathrm{N}\left( {0,{\sigma }^{2}}\right)\) 分布。如果得到的位置是允许的则接受结果,否则留在原地不动。

\section*{3.7.3 Gibbs 抽样}

一般的 MH 抽样每一步首先进行尝试运动, 然后根据新的状态是否靠近目标分布来接受或拒绝试抽样点, 所以可能会存在多次的无效尝试, 效率较低。

Gibbs 抽样是另外一种 MCMC 方法, 它仅在坐标轴方向尝试转移, 用当前点的条件分布决定下一步的试抽样分布, 所有试抽样都被接受, 不需要拒绝, 所以效率可以更高。

设状态用 \(\mathbf{x} = \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 表示,设目标分布为 \(\pi \left( \mathbf{x}\right)\) ,用 \({\mathbf{x}}_{\left( -i\right) }\) 表示 \(\left( {{x}_{1},\ldots ,{x}_{i - 1}}\right.\) , \(\left. {{x}_{i + 1},\ldots ,{x}_{n}}\right)\) ,假设 \(\pi \left( \cdot \right)\) 的条件分布 \(p\left( {{x}_{i} \mid  {\mathbf{x}}_{\left( -i\right) }}\right)\) 都能够比较容易地抽样。

Gibbs 抽样每一步从条件分布中抽样, 可以轮流从每一分量抽样, 这样的算法称为系统扫描 Gibbs 抽样算法:

从 \(\pi \left( \mathbf{x}\right)\) 的取值区域中任意取一个初值 \({\mathbf{X}}^{\left( 0\right) }\)

for \(\left( {t\text{ in }0 : \left( {N - 1}\right) }\right) \{\)

for \(\left( {i\text{ in }1 : n}\right) \{\)

从条件分布 \(p\left( {{x}_{i} \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 中抽取 \({X}_{i}^{ * }\)

\}

令 \({\mathbf{X}}^{\left( t + 1\right) } \leftarrow  \left( {{X}_{1}^{ * },\ldots ,{X}_{n}^{ * }}\right)\)

\}

从条件分布抽样的次序也可以是随机选取各个分量,这样的算法称为随机扫描 Gibbs 抽样算法:

从 \(\pi \left( \mathbf{x}\right)\) 的取值区域中任意取一个初值 \({\mathbf{X}}^{\left( 0\right) }\)

for \(\left( {t\text{ in }0 : \left( {N - 1}\right) }\right) \{\)

按概率 \(\mathbf{\alpha } = \left( {{\alpha }_{1},\ldots ,{\alpha }_{n}}\right)\) 从 \(\left( {1,\ldots ,n}\right)\) 中随机抽取下标 \(i\)

从条件分布 \(p\left( {{x}_{i} \mid  {\mathbf{X}}_{\left( -i\right) }^{\left( t\right) }}\right)\) 中抽取 \({X}_{i}^{ * }\) ,令 \({\mathbf{X}}^{\left( t + 1\right) } = \left( {{X}_{1}^{\left( t\right) },\ldots ,{X}_{i - 1}^{\left( t\right) },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\)

\}

其中下标的抽样概率 \(\alpha\) 为事先给定。

容易看出,无论采用系统扫描还是随机扫描的 Gibbs 抽样,如果 \({\mathbf{X}}^{\left( t\right) }\) 服从目标分布,则 \({\mathbf{X}}^{\left( t + 1\right) }\) 也服从目标分布。以系统扫描方法为例,设在第 \(t + 1\) 步已经抽取了 \({X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * }\) ,令 \(\mathbf{Y} = \left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) ,设 \(\mathbf{Y} \sim  \pi \left( \cdot \right)\) 。下一步从 \(\pi \left( \cdot \right)\) 的边缘密度 \(p\left( {{x}_{i} \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 抽取 \({X}_{i}^{ * }\) ,则 \({\mathbf{Y}}^{ * }\overset{\bigtriangleup }{ = }\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)\) 的分布密度在 \({\mathbf{Y}}^{ * }\) 处的值为

\[
p\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

\[
= p\left( {{X}_{i}^{ * } \mid  {X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right) p\left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

\[
= \pi \left( {{X}_{1}^{ * },\ldots ,{X}_{i - 1}^{ * },{X}_{i}^{ * },{X}_{i + 1}^{\left( t\right) },\ldots ,{X}_{n}^{\left( t\right) }}\right)
\]

即 \(\mathbf{Y} \sim  \pi \left( \cdot \right)\) 则 \({\mathbf{Y}}^{ * } \sim  \pi \left( \cdot \right)\) 。

例 3.7.4. 设目标分布为二元正态分布,设 \(\mathbf{X} \sim  \pi \left( \mathbf{x}\right)\) 为

\[
\mathrm{N}\left\{  {\left( \begin{array}{l} 0 \\  0 \end{array}\right) ,\left( \begin{array}{ll} 1 & \rho \\  \rho & 1 \end{array}\right) }\right\}
\]

采用系统扫描 Gibbs 抽样方案, 每一步的迭代为,

\[
\text{ 抽取 }{X}_{1}^{\left( t + 1\right) } \mid  {X}_{2}^{\left( t\right) } \sim  \mathrm{N}\left( {\rho {X}_{2}^{\left( t\right) },1 - {\rho }^{2}}\right)
\]

\[
\text{ 抽取 }{X}_{2}^{\left( t + 1\right) } \mid  {X}_{1}^{\left( t + 1\right) } \sim  \mathrm{N}\left( {\rho {X}_{1}^{\left( t + 1\right) },1 - {\rho }^{2}}\right)
\]

递推可得

\[
\left( \begin{array}{l} {X}_{1}^{\left( t\right) } \\  {X}_{2}^{\left( t\right) } \end{array}\right)  \sim  \mathrm{N}\left\{  {\left( \begin{matrix} {\rho }^{{2t} - 1}{X}_{2}^{\left( 0\right) } \\  {\rho }^{2t}{X}_{2}^{\left( 0\right) } \end{matrix}\right) ,\left( \begin{matrix} 1 - {\rho }^{{4t} - 2} & \rho  - {\rho }^{{4t} - 1} \\  \rho  - {\rho }^{{4t} - 1} & 1 - {\rho }^{4t} \end{matrix}\right) }\right\}   \tag{3.96}
\]

当 \(t \rightarrow  \infty\) 时, \(\left( {{X}_{1}^{\left( t\right) },{X}_{2}^{\left( t\right) }}\right)\) 的期望与目标分布期望之差为 \(O\left( {\left| \rho \right| }^{2t}\right)\) ,方差与目标分布方差之差为 \(O\left( {\left| \rho \right| }^{4t}\right)\) 。 例 3.7.5. 设目标分布为

\[
\pi \left( {x,y}\right)  \propto  \left( \begin{array}{l} n \\  x \end{array}\right) {y}^{x + \alpha  - 1}{\left( 1 - y\right) }^{n - x + \beta  - 1},x = 0,1,\ldots ,n,0 \leq  y \leq  1, \tag{3.97}
\]

则 \(X\left| {Y \sim  \mathrm{B}\left( {n,y}\right) ,Y}\right| X \sim  \operatorname{Beta}\left( {x + \alpha ,n - x + \beta }\right)\) 。易见 \(Y\) 的边缘分布为 \(\operatorname{Beta}\left( {\alpha ,\beta }\right)\) 。可以用 Gibbs 抽样方法模拟生成(X, Y)的样本链。

例 3.7.6. 在 Gibbs 抽样中, 每次变化的可以不是单个的分量, 而是两个或多个分量。例如, 设某个试验有 \(r\) 种不同结果,相应概率为 \(\mathbf{p} = \left( {{p}_{1},\ldots ,{p}_{r}}\right)\) (其中 \(\left. {\mathop{\sum }\limits_{{i = 1}}^{r}{p}_{i} = 1}\right)\) ,独立重复试验 \(n\) 次,各个结果出现的次数 \(\mathbf{X} = \left( {{X}_{1},\ldots ,{X}_{r}}\right)\) 服从多项分布。设 \(A = \left\{  {{X}_{1} \geq  1,\ldots ,{X}_{r} \geq  1}\right\}\) , 假设 \(P\left( A\right)\) 概率很小,要在条件 \(A\) 下对 \(\mathbf{X}\) 抽样,如果先生成 \(\mathbf{X}\) 的无条件样本再舍弃不符合条件 \(A\) 的部分则效率太低,可以采用如下的 Gibbs 抽样方法。

首先,任取初值 \({\mathbf{X}}^{\left( 0\right) }\) ,如 \({\mathbf{X}}^{\left( 0\right) } = \left( {1,\ldots ,1}\right)\) 。假设已生成了 \({\mathbf{X}}^{\left( t\right) }\) ,下一步首先从 \(\left( {1,\ldots ,r}\right)\) 中随机地抽取两个下标(i, j),令 \(s = {X}_{i}^{\left( t\right) } + {X}_{j}^{\left( t\right) }\) ,在给定 \({X}_{k},k \neq  i,j\) 为 \({\mathbf{X}}^{\left( t\right) }\) 对应元素的条件下, \(\left( {{X}_{i},{X}_{j}}\right)\) 的条件分布实际是 \(\left( {{X}_{i},{X}_{j}}\right)\) 在 \({X}_{i} + {X}_{j} = s\) 以及 \({X}_{i} \geq  1,{X}_{j} \geq  1\) 条件下的分布。于是,在以上条件下, \({X}_{i}\) 服从 \(\mathrm{B}\left( {s,{p}_{i}/\left( {{p}_{i} + {p}_{j}}\right) }\right)\) 分布限制在 \(1 \leq  {X}_{i} \leq  s - 1\) 条件下的分布,即

\[
{q}_{k}\overset{\bigtriangleup }{ = }P\left( {{X}_{i} = k \mid  {X}_{i} + {X}_{j} = s,{X}_{i} \geq  1,{X}_{j} \geq  1}\right)
\]

\[
= \frac{{C}_{s}^{k}{\left( \frac{{p}_{i}}{{p}_{i} + {p}_{j}}\right) }^{k}{\left( \frac{{p}_{j}}{{p}_{i} + {p}_{j}}\right) }^{s - k}}{1 - {\left( \frac{{p}_{j}}{{p}_{i} + {p}_{j}}\right) }^{s} - {\left( \frac{{p}_{i}}{{p}_{i} + {p}_{j}}\right) }^{s}},k = 1,2,\ldots ,s - 1.
\]

要生成这样的 \({X}_{i}\) 的抽样只要用生成离散型随机数的逆变换法。设抽取的 \({X}_{i}\) 值为 \({X}_{i}^{ * }\) ,取 \(\left( {{X}_{i}^{\left( t + 1\right) },{X}_{j}^{\left( t + 1\right) }}\right)  = \left( {{X}_{i}^{ * },s - {X}_{i}^{ * }}\right)\) ,取 \({\mathbf{X}}^{\left( t + 1\right) }\) 的其它元素为 \({\mathbf{X}}^{\left( t\right) }\) 的对应元素。如此重复就可以生成所需的 \(\mathbf{X}\) 在条件 \(A\) 下的抽样链。

\section*{3.7.4 MCMC 计算软件 *}

MCMC 是贝叶斯统计计算中最常用的计算工具。OpenBUGS 是一个成熟的 MCMC 计算开源软件 (见 Lunn et al(2009) \({}^{\left\lbrack  {29}\right\rbrack  }\) , Cowles(2013) \({}^{\left\lbrack  {16}\right\rbrack  }\) ,另一个类似的有关软件是 WinBUGS \({}^{\left\lbrack  {30}\right\rbrack  }\) ), 能够进行十分复杂的贝叶斯模型的计算,可以在 \(\mathrm{R}\) 中直接调用 OpenBUGS 进行计算。

OpenBUGS 采用 Gibbs 抽样方法从贝叶斯后验分布中抽样, 用户只需要指定先验分布和似然函数以及观测数据、已知参数, 以及并行地生成多少个马氏链、链的一些初值、运行步数。软件自动计算 Gibbs 抽样所需的条件分布, 产生马氏链, 并可以用图形和数值辅助判断收敛性, 给出后验推断的概括统计。

在 \(\mathrm{R}\) 中通过 BRugs 包调用 OpenBUGS 的功能。BRugs 用三类输入文件指定一个贝叶斯模型, 第一类文件指定似然函数和参数先验密度, 第二类文件指定已知参数、样本值, 第三类文件指定马氏链初值(并行产生多个链时需要指定多组初值)。 \(\mathrm{R}\) 的 coda 软件包可以帮助对 MCMC 抽样结果进行分析和诊断。下面用一个简单例子介绍在 \(\mathrm{R}\) 中用 BRugs 和 OpenBUGS 从贝叶斯后验中抽样的基本步骤。

例 3.7.7. 对例3.7.2,用 \(\mathrm{R}\) 的 BRugs 包调用 OpenBUGS 来计算。设 \(\left( {{X}_{1},\ldots ,{X}_{5}}\right)\) 的观测值为(74,85,69,17,5)。

OpenBUGS 用模型文件描述随机变量分布和对参数的依赖关系, 以及参数之间的关系。 首先建立如下模型文件,保存在文件 pfl-model.txt 中:

\HRule

model

\{

\hspace*{1em} \(p\left\lbrack  1\right\rbrack   <  - 1/3\)

\hspace*{1em} \(p\left\lbrack  2\right\rbrack   <  - \left( {1 - b}\right) /3\)

\hspace*{1em} \(p\left\lbrack  3\right\rbrack   <  - \left( {1 - 2 * b}\right) /3\)

\hspace*{1em} \(p\left\lbrack  4\right\rbrack   <  - 2 * b/3\)

\hspace*{1em} \(p\left\lbrack  5\right\rbrack   <  - b/3\)

\hspace*{1em} b <- b2 / 2

\hspace*{1em} b2 \~ dbeta(1,1)

\hspace*{1em} x[1:5] \~ dmulti(p[1:5], N)

\}

\HRule

文件中用向左的箭头表示确定性的关系,用方括号加序号表示下标,用 \(\sim\) 表示左边的变量服从右边的分布。这里, \(b\) 为参数 \(\beta ,\mathrm{b}2 = {2b}\) 服从 \(\operatorname{Beta}\left( {1,1}\right)\) 分布即 \(\mathrm{U}\left( {0,1}\right)\) 分布,于是 \(\beta\) 有先验分布 \(\mathrm{U}\left( {0,{0.5}}\right)\) 。 \(\mathrm{x}\left\lbrack  {1 : 5}\right\rbrack\) 表示向量 \(\left( {{x}_{1},\ldots ,{x}_{5}}\right)\) ,服从多项分布,参数为 \(\left( {{p}_{1},\ldots ,{p}_{5}}\right)\) ,试验次数为 \(N\) 。 \(N\) 在模型文件中没有指定,将在数据文件中给出。

建立如下的数据文件, 保存在文件 pfl-data.txt 中:

\HRule

list(x=c(74, 85, 69, 17, 5),

\hspace*{1em} N=250)

\HRule

这样的数据文件是 \(\mathrm{R}\) 软件的列表格式,列表中的标量和向量为通常的 \(\mathrm{R}\) 程序写法,矩阵用如

\HRule

list(M=structure(.Data=c(1,2,3,4,5,6),.Dim=c(3,2)))

\HRule

表示,其中. Dim 给出矩阵的行、列数 \(\left( {3 \times  2}\right)\) ,. Data 给出按行排列的所有元素(第一行为 1,2,第二行为3,4,第三行为5,6)。

OpenBUGS 需要用户指定各个链的要抽样的参数的初值, 这里我们要抽样的参数是 b, 但 \(\mathrm{b}\) 是由 \(\mathrm{b}2\) 计算得到的,所以对 \(\mathrm{b}2\) 设置初值。设有如下两个初值文件,不同链的初值应

尽可能不同, 文件 pfl-inits1.txt 内容为:

\HRule

list(b2=0.2)

\HRule

文件 pfl-inits2.txt 内容为:

\HRule

list(b2=0.8)

\HRule

在 \(\mathrm{R}\) 中,首先调用 BRugs 软件包:

\HRule

require (BRugs)

\HRule

然后, 读入并检查模型文件:

\HRule

modelCheck('pfl-model.txt')

\HRule

读入数据文件:

\HRule

modelData ( ' pfl-data.txt ' )

\HRule

下面, 对模型和数据进行编译, 得到抽样方案, 下面的语句要求并行运行两个链:

\HRule

modelCompile(numChains=2)

\HRule

准备迭代地生成 MCMC 抽样了, 首先指定初值:

\HRule

modelInits(c('pfl-inits1.txt', 'pfl-inits2.txt'))

\HRule

下面, 先试验性地运行 1000 次, 作为老化期:

\HRule

modelUpdate (1000)

\HRule

现在才指定抽样要输出那些随机变量的随机数:

\HRule

samplesSet \(\left( {c\left( {{}^{\prime }{b}^{\prime },{}^{\prime }p\left\lbrack  {1 : 5}\right\rbrack  {}^{\prime }}\right) }\right)\)

\HRule

现在可以抽样了, 指定运行 10000 次, 两个链并行运行:

\HRule

modelUpdate (10000)

\HRule

得到抽样后,可以把抽样的结果保存在 \(\mathrm{R}\) 的变量中,比如

\HRule

b2chains <- samplesHistory ( 'b', plot=FALSE )

\HRule

得到一个列表,有唯一的元素 \(\mathrm{b}\) ,为 \(2 \times  {10000}\) 的矩阵,每行是一个链的记录。sampleHistory 可以抽样链的曲线图, 如

samplesHistory('b')

OpenBUGS 提供了一系列的简单统计和收敛诊断功能。如下程序列出各抽样变量的简单统计:

\HRule

print ( samplesStats ( "*" ) )

\HRule

结果为:

\HRule

\hspace*{1em} <table><tr><td/><td>mean</td><td>sd</td><td>MC\_error</td><td>val2.5pc</td><td>median</td><td>val97.5pc</td><td>start</td><td>sample</td></tr><tr><td>b</td><td>0.08757</td><td>0.016900</td><td>1.168e-04</td><td>0.05737</td><td>0.08668</td><td>0.12330</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  2\right\rbrack\) </td><td>0.30410</td><td>0.005632</td><td>3.892e-05</td><td>0.29230</td><td>0.30440</td><td>0.31420</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  3\right\rbrack\) </td><td>0.27500</td><td>0.011260</td><td>7.784e-05</td><td>0.25120</td><td>0.27550</td><td>0.29510</td><td>1001</td><td>20000</td></tr><tr><td>p [4]</td><td>0.05838</td><td>0.011260</td><td>7.784e-05</td><td>0.03824</td><td>0.05778</td><td>0.08217</td><td>1001</td><td>20000</td></tr><tr><td> \(\mathrm{p}\left\lbrack  5\right\rbrack\) </td><td>0.02919</td><td>0.005632</td><td>3.892e-05</td><td>0.01912</td><td>0.02889</td><td>0.04109</td><td>1001</td><td>20000</td></tr></table>

\HRule

统计使用所有链的数据。可以看出, \(\beta\) 的后验均值为 0.08757。表中的 \({\mathrm{{MC}}}_{ - }\) error 表示估计后验均值时由于随机模拟导致的误差的标准差的估计, 这个标准差估计针对抽样自相关性进行了校正。在老化期之后的运行次数越多 MC\_error 越小, 一个常用的经验规则是保证 MC\_error 小于后验标准差的 \(5\%\) ,这里的结果提示还需要更多的运行次数。val \({2.5}\mathrm{{pc}}\) 和 val \({97.5}\mathrm{{pc}}\) 是抽样的后验分布的 2.5\% 和 97.5\% 分位数的估计值,由此得到 \(\beta\) 的水平 95\% 的可信区间 (credible interval) 为(0.05734,0.1233)。

如下程序画出 \(\mathrm{b}\) 的后验密度估计:

\HRule

samplesDensity ( ' b ', mfrow=c ( 1 , 1 ) )

\HRule

如下程序画抽样的 \(\mathrm{b}\) 的自相关函数估计:

\HRule

samplesAutoC('b', 1, mfrow=c(1,1))

\HRule

当自相关函数很大而且衰减缓慢时生成的抽样链的效率较低。本例中的自相关函数基本表现为不相关列。

为了检查链是否收敛, OpenBUGS 提供了 BGR 统计量图:

\HRule

samplesBgr ( ' b ', mfrow=c (1,1 ) )

\HRule

BGR 统计量的原理是考虑并行运行的每个链内部的变化情况, 以及把所有的链合并在一起的变化情况, 对这两种变化情况进行比较, 当链收敛时, 每个链内部的变化情况应该与合并在一起的变化情况很类似。类似于单因素方差分析中组间平方和与组内平方和的比较。当 BGR 图中的红色线接近于 1 并且三条线都保持稳定时就提示链收敛了。

\section*{3.8 序贯重要抽样 *}

MCMC 是目前广泛使用的随机模拟方法, 其中 Gibbs 抽样方法 (见 §3.7.3) 在确定了各个分量的条件分布后可以轮流产生各个分量的抽样。序贯重要抽样方法是重要抽样法(见 §3.2.4)的一种推广,其做法与 Gibbs 抽样方法有些相似,也是每次从一个分量抽样。

回顾多维随机变量抽样的条件分布法 (见 §2.3.1)。设随机向量 \(\mathbf{X}\) 的密度 \(\pi \left( \mathbf{x}\right)  =\)  \(\pi \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)\) 可以逐次地分解为条件密度乘积

\[
\pi \left( {{x}_{1},{x}_{2},\ldots ,{x}_{n}}\right)  = \pi \left( {x}_{1}\right) \pi \left( {{x}_{2} \mid  {x}_{1}}\right) \pi \left( {{x}_{3} \mid  {x}_{1},{x}_{2}}\right) \cdots \pi \left( {{x}_{n} \mid  {x}_{1},\ldots ,{x}_{n - 1}}\right)
\]

则可以用条件分布法产生 \(\mathbf{X}\) 的抽样。当数据是时间序列或者可以依次增加对 \(\pi\) 的信息(比如 \(\pi\) 是 Bayes 后验分布) 时这种方法很自然。但是,条件分布 \(\pi \left( {{x}_{t} \mid  {x}_{1},\ldots ,{x}_{t - 1}}\right)\) 可能是难以得到的或难以抽样的。为此,采用重要抽样思想: 取一系列辅助分布 \({\pi }_{t}\left( {\mathbf{x}}_{t}\right)  = {\pi }_{t}\left( {{x}_{1},\ldots ,{x}_{t}}\right) ,t =\)  \(1,\ldots ,n\) ,使其近似于 \(\left( {{X}_{1},\ldots ,{X}_{t}}\right)\) 的分布, \({\pi }_{n} = \pi\) ,各 \({\pi }_{t}\) 可以分别有一个未知的归一化常数。对 \(t = 1,2,\ldots ,n\) 用重要抽样法逐次抽取 \({X}_{1},{X}_{2},\ldots ,{X}_{n}\) ,使得 \(\left( {{X}_{1},\ldots ,{X}_{t}}\right)\) 是关于 \({\pi }_{t}\) 适当加权的样本。抽取 \({X}_{t}\) 时一般是给定 \({X}_{1},\ldots ,{X}_{t - 1}\) 后从一个试抽样分布 \({g}_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)\) 抽取。适当计算权重就可以得到关于 \(\pi\) 适当加权的抽样 \({\mathbf{X}}_{n} = \left( {{X}_{1},{X}_{2},\ldots ,{X}_{n}}\right)\) 。

这种抽样方法叫做序贯重要抽样 (sequential importance sampling, SIS), 算法如下:

置 \(t \leftarrow  1\) ,从 \({g}_{1}\left( \cdot \right)\) 抽取 \({X}_{1}\) ,置 \({W}_{1} \leftarrow  {\pi }_{1}\left( {X}_{1}\right) /{g}_{1}\left( {X}_{1}\right)\) 。

for \(\left( {t\text{ in }2 : n}\right) \{\)

从 \({g}_{t}\left( {{x}_{t} \mid  {X}_{1},\ldots ,{X}_{t - 1}}\right)\) 抽取 \({X}_{t}\) ,记 \({\mathbf{X}}_{t} = \left( {{X}_{1},\ldots ,{X}_{t - 1},{X}_{t}}\right)\) ;

计算步进权重

\[
{U}_{t} \leftarrow  \frac{{\pi }_{t}\left( {\mathbf{X}}_{t}\right) }{{\pi }_{t - 1}\left( {\mathbf{X}}_{t - 1}\right) {g}_{t}\left( {{X}_{t} \mid  {\mathbf{X}}_{t - 1}}\right) }
\]

令 \({W}_{t} \leftarrow  {W}_{t - 1}{U}_{t}\) ;

\}

输出 \(\left( {{\mathbf{X}}_{n},{W}_{n}}\right)\) 为关于 \(\pi \left( \mathbf{x}\right)\) 适当加权的样本。

SIS 一般同时独立进行 \(N\) 组,得到 \({\left\{  \left( {\mathbf{X}}_{n}^{\left( j\right) },{W}_{n}^{\left( j\right) }\right) \right\}  }_{j = 1}^{N}\) ,每一组称为一个 “流” 或一个 “粒子”。 按照 SIS 步骤, 有

\[
{X}_{1} \sim  {g}_{1}\left( {x}_{1}\right) ,\;{w}_{1} = \frac{{\pi }_{1}\left( {X}_{1}\right) }{{g}_{1}\left( {X}_{1}\right) }
\]

\[
{X}_{2} \sim  {g}_{2}\left( {{x}_{2} \mid  {x}_{1}}\right) ,\;{U}_{2} = \frac{{\pi }_{2}\left( {\mathbf{X}}_{2}\right) }{{\pi }_{1}\left( {X}_{1}\right) {g}_{2}\left( {{X}_{2} \mid  {X}_{1}}\right) },
\]

\[
{W}_{2} = {W}_{1}{U}_{2} = \frac{{\pi }_{2}\left( {\mathbf{X}}_{2}\right) }{{g}_{1}\left( {X}_{1}\right) {g}_{2}\left( {{X}_{2} \mid  {X}_{1}}\right) }
\]

\[
{X}_{3} \sim  {g}_{3}\left( {{X}_{3} \mid  {\mathbf{X}}_{2}}\right) ,\;{U}_{3} = \frac{{\pi }_{3}\left( {\mathbf{X}}_{3}\right) }{{\pi }_{2}\left( {\mathbf{X}}_{2}\right) {g}_{3}\left( {{X}_{3} \mid  {\mathbf{X}}_{2}}\right) },
\]

\[
{W}_{3} = \frac{{\pi }_{3}\left( {\mathbf{X}}_{3}\right) }{{g}_{1}\left( {X}_{1}\right) {g}_{2}\left( {{X}_{2} \mid  {X}_{1}}\right) {g}_{3}\left( {{X}_{3} \mid  {\mathbf{X}}_{2}}\right) }
\]

......

\[
{X}_{n} \sim  {g}_{n}\left( {{x}_{n} \mid  {\mathbf{X}}_{n - 1}}\right) ,\;{W}_{n} = \frac{{\pi }_{n}\left( {\mathbf{X}}_{n}\right) }{{g}_{1}\left( {X}_{1}\right) {g}_{2}\left( {{X}_{2} \mid  {X}_{1}}\right) \cdots {g}_{n}\left( {{X}_{n} \mid  {\mathbf{X}}_{n - 1}}\right) }.
\]

记

\[
{g}_{t}\left( {\mathbf{x}}_{t}\right)  = {g}_{1}\left( {x}_{1}\right) {g}_{2}\left( {{x}_{2} \mid  {x}_{1}}\right) \cdots {g}_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)
\]

则 \({\mathbf{X}}_{t} \sim  {g}_{t}\left( \cdot \right)\) 而

\[
{W}_{t} = \frac{{\pi }_{t}\left( {\mathbf{X}}_{t}\right) }{{g}_{t}\left( {\mathbf{X}}_{t}\right) },
\]

因此 \(\left( {{\mathbf{X}}_{t},{W}_{t}}\right)\) 关于 \({\pi }_{t}\) 适当加权 \(\left( {t = 1,\ldots ,n}\right)\) 。注意 \({\pi }_{n} = \pi\) 故这样得到的 \(\left( {{\mathbf{X}}_{n},{W}_{n}}\right)\) 关于 \(\pi \left( \cdot \right)\) 适当加权。

试抽样分布的一种常见取法为 \({g}_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)  = {\pi }_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)\) 。

\section*{3.8.1 非线性滤波平滑}

SIS 方法可以用在很多统计模型的计算中,作为示例,考虑如下的非线性滤波平滑问题。

设不可观测的 “状态” \({X}_{t}\) 服从如下状态方程

\[
{X}_{t} \sim  {q}_{t}\left( {\cdot  \mid  {X}_{t - 1},\theta }\right) ,t = 1,2,\ldots ,n
\]

\({X}_{t}\) 的信息可以反映在可观测的 \({Y}_{t}\) 中,其关系满足如下观测方程

\[
{Y}_{t} \sim  {f}_{t}\left( {\cdot  \mid  {X}_{t},\phi }\right) ,t = 1,2,\ldots ,n
\]

已知 \({\mathbf{Y}}_{n} = \left( {{Y}_{1},\ldots ,{Y}_{n}}\right)  = {\mathbf{y}}_{n}\) 和 \(\theta ,\phi\) 后估计 \(\mathbf{X} = \left( {{X}_{1},\ldots ,{X}_{n}}\right)\) 的问题称为滤波平滑问题,只需要求后验分布 \(\pi \left( {\mathbf{x}}_{n}\right)  = {p}_{{\mathbf{X}}_{n} \mid  {\mathbf{Y}}_{n}}\left( {{\mathbf{x}}_{n} \mid  {\mathbf{y}}_{n}}\right)\) 。如果能从 \(\pi\) 大量抽样 \({\mathbf{X}}^{\left( j\right) },j = 1,\ldots ,N\) 则可以用随机模拟方法对 \({\mathbf{X}}_{n} = \left( {{X}_{1},\ldots ,{X}_{n}}\right)\) 的后验分布进行推断。下面用 SIS 方法产生 \(\pi\) 的样本。

记

\[
{\pi }_{t}\left( {\mathbf{x}}_{t}\right)  = {p}_{{\mathbf{X}}_{t} \mid  {\mathbf{Y}}_{t}}\left( {{\mathbf{x}}_{t} \mid  {\mathbf{y}}_{t}}\right) ,t = 1,2,\ldots ,n
\]

注意

\[
{\pi }_{t}\left( {\mathbf{x}}_{t}\right)  \propto  {f}_{t}\left( {{y}_{t} \mid  {x}_{t}}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right) {\pi }_{t - 1}\left( {\mathbf{x}}_{t - 1}\right)  \tag{3.98}
\]

取试抽样分布 \({g}_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)  = {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right)\) ,对 \(t = 1,{\pi }_{1}\left( {x}_{1}\right)  = {p}_{{X}_{1} \mid  {Y}_{1}}\left( {{x}_{1} \mid  {y}_{1}}\right)  \propto  {f}_{1}\left( {{y}_{1} \mid  {x}_{1}}\right) p\left( {x}_{1}\right)\) ,其中 \(p\left( {x}_{1}\right)\) 为 \({X}_{1}\) 的分布密度,抽取 \({X}_{1} \sim  {g}_{1}\left( {x}_{1}\right)\) ,如果可能应取 \({g}_{1}\left( {x}_{1}\right)  = {\pi }_{1}\left( {x}_{1}\right)\) 。

产生关于 \(\pi\) 适当加权的 \({\mathbf{X}}_{n}\) 抽样的 SIS 步骤如下:

置 \(t \leftarrow  1\) ,从 \({g}_{1}\left( {x}_{1}\right)\) 抽取 \({X}_{1}\) ,置 \({W}_{1} \leftarrow  {\pi }_{1}\left( {X}_{1}\right) /{g}_{1}\left( {X}_{1}\right)\) 。

for \(\left( {t\text{ in }2 : n}\right) \{\)

从 \({q}_{t}\left( {{x}_{t} \mid  {X}_{t - 1}}\right)\) 抽取 \({X}_{t}\) ,记 \({\mathbf{X}}_{t} = \left( {{X}_{1},\ldots ,{X}_{t - 1},{X}_{t}}\right)\) ;

令步进权重 \({U}_{t} \leftarrow  {f}_{t}\left( {{y}_{t} \mid  {X}_{t}}\right)\)

令 \({W}_{t} \leftarrow  {W}_{t - 1}{U}_{t}\) ;

\}

输出 \(\left( {{\mathbf{X}}_{n},{W}_{n}}\right)\) 为关于 \(\pi \left( \mathbf{x}\right)\) 适当加权的样本。

这种方法相当于用状态方程前进一步获得下一分量的抽样,用同时刻的观测值 \({y}_{t}\) 的似然作为步进权重。这样, \(t\) 步以后得到的 \({\mathbf{X}}_{t}\) 服从

\[
{g}_{t}\left( {\mathbf{x}}_{t}\right)  = {g}_{t - 1}\left( {\mathbf{x}}_{t - 1}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right)
\]

其中 \({g}_{1}\left( {x}_{1}\right)  \propto  {f}_{1}\left( {{y}_{1} \mid  {x}_{1}}\right)\) 。于是

\[
\frac{{\pi }_{t}\left( {\mathbf{x}}_{t}\right) }{{g}_{t}\left( {\mathbf{x}}_{t}\right) } = \frac{{f}_{t}\left( {{y}_{t} \mid  {x}_{t}}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right) {\pi }_{t - 1}\left( {\mathbf{x}}_{t - 1}\right) }{{g}_{t - 1}\left( {\mathbf{x}}_{t - 1}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right) }
\]

\[
= {f}_{t}\left( {{y}_{t} \mid  {x}_{t}}\right) \frac{{\pi }_{t - 1}\left( {\mathbf{x}}_{t - 1}\right) }{{g}_{t - 1}\left( {\mathbf{x}}_{t - 1}\right) },
\]

\[
{W}_{t} = {f}_{t}\left( {{y}_{t} \mid  {X}_{t}}\right) {W}_{t - 1} = \frac{{\pi }_{t}\left( {\mathbf{X}}_{t}\right) }{{g}_{t}\left( {\mathbf{X}}_{t}\right) },
\]

可见 \(\left( {{\mathbf{X}}_{t},{W}_{t}}\right)\) 关于 \({\pi }_{t}\) 适当加权, \(\left( {{\mathbf{X}}_{n},{W}_{n}}\right)\) 关于 \({\pi }_{n}\left( \mathbf{x}\right)  = {p}_{{\mathbf{X}}_{n} \mid  {\mathbf{Y}}_{n}}\left( {{\mathbf{x}}_{n} \mid  {\mathbf{y}}_{n}}\right)\) 适当加权。

设第 \(t\) 步的抽样为 \({\mathbf{X}}_{t}^{\left( i\right) },i = 1,\ldots ,N\) ,对应权重为 \({W}_{t}^{\left( i\right) },i = 1,\ldots ,n\) 。以上的方法在抽取 \(\left( {{X}_{1},\ldots ,{X}_{n}}\right)\) 时不考虑 \(\left( {{y}_{1},\ldots ,{y}_{n}}\right)\) 的具体取值,这样的试抽样分布虽然很容易抽样,但是效果很差,权重 \(\left\{  {{W}_{t}^{i},i = 1,\ldots ,N}\right\}\) 随着 \(t\) 的增加会把变得差异很大,以至于 \(N\) 个流中只有极少数流能起作用。

由(3.98)可见

\[
{\pi }_{t}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)  = {p}_{{X}_{t} \mid  {\mathbf{X}}_{t - 1},{Y}_{t}}\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)  \propto  {f}_{t}\left( {{y}_{t} \mid  {x}_{t}}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right) , \tag{3.99}
\]

如果能取试抽样分布 \({g}_{t}\left( {x}_{t}\right)  \propto  {f}_{t}\left( {{y}_{t} \mid  {x}_{t}}\right) {q}_{t}\left( {{x}_{t} \mid  {x}_{t - 1}}\right)\) 则每次抽取 \({X}_{t}\) 都利用了同期的观测值 \({y}_{t}\) 的信息,会大大改善抽样效率。更进一步,设 \({\pi }_{t + 1}\left( {\mathbf{x}}_{t + 1}\right)\) 关于 \({\mathbf{x}}_{t}\) 的边缘分布为 \({\pi }_{t,t + 1}\left( {\mathbf{x}}_{t}\right)\) ,如果在第 \(t\) 步能从 \({\pi }_{t,t + 1}\left( {\mathbf{x}}_{t}\right)\) 的条件分布 \(p\left( {{x}_{t} \mid  {\mathbf{x}}_{t - 1}}\right)\) 抽样,就可以利用 \({y}_{t},{y}_{t + 1}\) 的信息,抽样效率可以进一步改善。

另外一种改进的办法是再抽样,增加权重大的流,舍弃权重小的流。

\section*{3.8.2 再抽样}

如果试抽样分布选取不适当,最后的权重可能会差别很大,体现在权重 \(\left\{  {{W}_{t}^{\left( i\right) },i = }\right.\)  \(1,\ldots ,N\}\) 的样本变异系数很大,称为权重偏斜严重。这时,权重小的流基本不起作用。出现这样的情况时,可以把一些权重太小的流舍弃而增加权重大的流,这样的技术称为再抽样。

\section*{简单随机再抽样}

设进行 SIS 时在时刻 \(t\) 已经得到了 \(N\) 个部分的流 \(\left\{  {{\mathbf{X}}_{t}^{\left( i\right) },i = 1,\ldots ,N}\right\}\) 以及相应的权重 \(\left\{  {{W}_{t}^{\left( i\right) },i = 1,\ldots ,N}\right\}\) 。可以在每一步都按照权重对流再抽样,也可以仅当权重偏斜严重时才进行再抽样。如果在第 \(t\) 步再抽样,只要以正比于 \({W}_{t}^{\left( i\right) }\) 的概率从 \(\left\{  {{\mathbf{X}}_{t}^{\left( j\right) },j = 1,\ldots ,N}\right\}\) 中抽取 \({\mathbf{X}}_{t}^{\left( i\right) }\) ,独立有放回地抽取 \(N\) 个,记作 \(\left\{  {{\mathbf{X}}_{t}^{*\left( i\right) },i = 1,\ldots ,N}\right\}\) ,并把权重调整为相等的 \({W}_{t}^{*\left( i\right) } = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{W}_{t}^{\left( j\right) },i = 1,\ldots ,N\) 。这样的方法称为简单随机再抽样。这样再抽样后各个流不再是独立的。

\section*{剩余再抽样}

为了达到以上的简单随机抽样的效果, 还可以把大权重的的流直接复制多份, 小权重的流仅以一定比例保留,其它舍弃。这种方法称为剩余再抽样 (residual resampling),其计算量更小而且模拟误差更小。算法描述如下。如果在第 \(t\) 步需要再抽样,则计算 \({\bar{W}}_{t} = \frac{1}{N}\mathop{\sum }\limits_{{j = 1}}^{N}{W}_{t}^{\left( j\right) }\) , 然后按如下做法从 \(N\) 个流中重新抽取 \(N\) 个。首先,对 \(i = 1,\ldots ,N\) ,直接保留 \({k}_{i} = \left\lbrack  {{W}_{t}^{\left( i\right) }/{\bar{W}}_{t}}\right\rbrack\) 份 \({\mathbf{X}}_{t}^{\left( i\right) }\left( {\lbrack  \cdot  \rbrack \text{ 表示向下取整); 其次,令 }{N}_{r} = N - \mathop{\sum }\limits_{{i = 1}}^{N}{k}_{i}}\right)\) 为缺额个数,随机有放回地按照正比于 \(\frac{{W}_{t}^{\left( i\right) }}{{W}_{t}} - {k}_{i}\) (取整后的小数部分) 的概率从 \(\left\{  {{\mathbf{X}}_{t}^{\left( j\right) },j = 1,\ldots ,N}\right\}\) 中抽取 \({\mathbf{X}}_{t}^{\left( i\right) }\) ,共抽取 \({N}_{r}\) 个。这样得到了 \(N\) 个新的流,记作 \(\left\{  {{\mathbf{X}}_{t}^{*\left( i\right) },i = 1,\ldots ,N}\right\}\) ,并调整其权重为相等的 \({W}_{t}^{*\left( i\right) } = {\bar{W}}_{t}\) , \(i = 1,\ldots ,N\) 。这种方法的第一步把权重大的流直接复制了若干份,第二步对按剩余的权重再抽取到满 \(N\) 个流为止。

\section*{舍选控制再抽样}

简单随机再抽样和剩余再抽样都使得结果各个流不再独立。另外一种想法是采用 \(§{3.2.4}\) 中介绍的舍选控制方法, 对权重小的流从头重新抽样并适当调整权重。首先, 设定若干个要执行再抽样的时间点 \(0 < {t}_{1} < \cdots  < {t}_{k} \leq  n\) ,以及相应的权重阈值 \({c}_{1},\ldots ,{c}_{k}\) 。在 \(t = {t}_{j}\) 时,进行舍选控制再抽样。若流 \(i\) 的权重 \({W}_{t}^{\left( i\right) } \geq  {c}_{j}\) ,则保留此流 \({\mathbf{X}}_{t}^{\left( i\right) }\) 和权重 \({W}_{t}^{\left( i\right) }\) 不变; 若流 \(i\) 的权重 \({W}_{t}^{\left( i\right) } < {c}_{j}\) ,则以概率 \({W}_{t}^{\left( i\right) }/{c}_{j}\) 保留此流,如果决定保留,则修改其权重为 \({c}_{j}\) ,如果决定舍弃此流,则从 \(t = 1\) 重新生成这个流,同样也需要经过 \({t}_{1},\ldots ,{t}_{j}\) 处的舍选判断,如果被舍弃就还从 \(t = 1\) 重新开始,直到被接受。

\section*{部分舍选控制再抽样}

如果从头开始重新抽样的情况发生比较多则模拟的效率会比较低,为此可以采用如下的部分舍选控制再抽样的 SIS 方法。

首先,设定若干个要执行再抽样的时间点 \(0 < {t}_{1} < \cdots  < {t}_{k} \leq  n\) ,以及相应的权重阈值 \({c}_{1},\ldots ,{c}_{k}\) 。在 \(t = {t}_{j}\) 时,进行部分舍选控制再抽样。若流 \(i\) 的权重 \({W}_{t}^{\left( i\right) } \geq  {c}_{j}\) ,则保留此流 \({\mathbf{X}}_{t}^{\left( i\right) }\) 和权重 \({W}_{t}^{\left( i\right) }\) 不变。若流 \(i\) 的权重 \({W}_{t}^{\left( i\right) } < {c}_{j}\) ,则以概率 \({W}_{t}^{\left( i\right) }/{c}_{j}\) 保留此流,如果决定保留,则修改其权重为 \({c}_{j}\) 。如果决定舍弃此流,则不是从头重新生成这个流,而是按照概率正比于 \({W}_{{t}_{j - 1}}^{\left( s\right) }\) 从 \(\left\{  {{\mathbf{X}}_{{t}_{j - 1}}^{\left( s\right) },s = 1,\ldots ,N}\right\}\) 中随机抽取一个替换原来的 \({\mathbf{X}}_{{t}_{j - 1}}^{\left( i\right) }\) ,用 \({t}_{j - 1}\) 时的权重 \(\left\{  {W}_{{t}_{j - 1}}^{\left( s\right) }\right\}\) 的平均值 \({\bar{W}}_{{t}_{j - 1}}\) 代替原来的权重 \({W}_{{t}_{j - 1}}^{\left( i\right) }\) ,然后继续按照 SIS 标准步骤将此流经 \(t = {t}_{j - 1} + 1,\ldots ,{t}_{j}\) 延伸得到新的 \({\mathbf{X}}_{{t}_{j}}^{\left( i\right) }\) 和权重 \({W}_{{t}_{j}}^{\left( i\right) }\) ,然后再进行 \({t}_{j}\) 处的舍选控制,如果被舍弃就再从 \({t}_{j - 1}\) 处随机再抽样后继续,直到在 \({t}_{j}\) 处被接受。

关于序贯重要抽样更详细的讨论参见 \(\operatorname{Liu}{\left( {2001}\right) }^{\left\lbrack  {28}\right\rbrack  }\) 。

\section*{习题三}

1. 设 \({\sigma }_{j},j = 1,2,\ldots ,m\) 为 \(m\) 个正实数,

\[
f\left( {{\alpha }_{1},{\alpha }_{2},\ldots ,{\alpha }_{m}}\right)  = \frac{{\sigma }_{1}^{2}}{{\alpha }_{1}} + \frac{{\sigma }_{2}^{2}}{{\alpha }_{2}} + \cdots  + \frac{{\sigma }_{m}^{2}}{{\alpha }_{m}},\left( {{\alpha }_{1},\ldots ,{\alpha }_{m}}\right)  \in  {\left( 0,1\right) }^{m},
\]

则在 \({\alpha }_{1} + {\alpha }_{2} + \cdots  + {\alpha }_{m} = 1\) 条件下 \(f\left( {{\alpha }_{1},{\alpha }_{2},\ldots ,{\alpha }_{m}}\right)\) 的最小值点为

\[
{\alpha }_{j} = \frac{{\sigma }_{j}}{\mathop{\sum }\limits_{{k = 1}}^{m}{\sigma }_{k}},j = 1,2,\ldots ,m
\]

最小值为 \({\left( {\sigma }_{1} + \cdots  + {\sigma }_{m}\right) }^{2}\) 。

2. 考虑定积分

\[
I = {\int }_{-1}^{1}{e}^{x}{dx} = e - {e}^{-1}.
\]

(1)用随机模拟方法计算定积分 \(I\) ,分别用随机投点法、平均值法、重要抽样法和分层抽样法计算。

(2)设估计结果为 \(\widehat{I}\) ,如果需要以 \({95}\%\) 置信度保证计算结果精度在小数点后三位小数, 这四种方法分别需要计算多少次被积函数值?

( 3 )用不同的随机数种子重复以上的估计 \(B\) 次,得到 \({\widehat{I}}_{j},j = 1,2,\ldots ,B\) ,由此估计 \(\widehat{I}\) 的抽样分布方差,与(2)的结果进行验证。

(4) 称

\[
\operatorname{MAE}\left( \widehat{I}\right)  = E\left| {\widehat{I} - I}\right|
\]

为 \(\widehat{I}\) 的平均绝对误差。从 (3) 得到的 \({\widehat{I}}_{j},j = 1,2,\ldots ,B\) 中估计 \(\operatorname{MAE}\left( \widehat{I}\right)\) 。比较这四种积分方法的平均绝对误差大小。
\end{document}